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ABSTRACT 


The  estimation  of  time  varying  spectra  is  a  complicated  one.  The  use  of  classical 
techniques  coupled  with  the  local  stationarity  assumption  is  met  with  only  moderate 
success.  Of  the  many  time-frequency  distribution  functions  used  in  the  signal  analysis, 
none  present  fully  satisfactory  spectra.  The  performance  of  the  spectrogram,  Instanta¬ 
neous  Power  Spectra  (IPS)  the  Wigner-Ville  distribution  (WD)  and  various  aspects  of 
the  Rihaczek  distribution  (RD)  for  a  variety  of  signal  nonstationarities  are  compared. 
WD  has  the  most  narrow  main-lobes  but  suffers  from  spectral  cross-terms.  IPS,  the  real 
part  of  the  RD  consistently  shows  a  bro  jened  main-lobe  without  cross-terms.  The 
squared  magnitude  of  the  RD  places  sharp  peaks  along  the  crest  of  the  main-lobe  and 
is  otherwise  very  similar  to  IPS.  The  imaginary  part  of  the  RD  shows  a  sensitivity  to 
discontinuous  frequency  changes  i.e.,  frequencv  shut '  eying. 
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I.  INTRODUCTION 


The  analysis  of  stationary  spectra  is  a  \v  ill-defined  problem.  Although  from  a  the¬ 
oretical  point  of  view  a  true  signal  spectrum  can  only  be  defined  in  terms  of  infinite  du¬ 
ration  data,  spectral  estimates  resulting  from  the  analysis  of  finite  duration  data  have 
proven  very  useful.  In  the  classical  estimation  problem,  one  assumes  that 

•  the  data  length  is  finite,  and 

•  the  random  process  from  which  it  originates  is  at  least  wide-sense  stationary. 

Armed  with  these  assumptions  the  behavior  of  spectra  derived  in  this  manner  can  be 
accurately  predicted.  The  distortion  incurred  when  analyzing  a  finite  amount  of  data  can 
be  kept  at  a  minimum  given  the  specific  estimation  problem.  But  what  happens  if  the 
signal  is  not  at  least  wide  sense  stationary?  Signals  commonly  encountered  in  the  real 
world  are  not  stationary;  they  vary  in  time.  Either  in  amplitude  or  frequency  content 
or  possibly  both,  experimental  data  is  rarely  truly  stationary.  To  better  describe  the 
variable  random  process,  a  time  dependence  must  be  included.  All  the  theoretical  results 
and  even  the  practical  application  to  finite  duration  data  assume  that  ergodicity  applies. 
Clearly  a  nonstationarv  signal  is  not  an  ergodic  one,  i.e.,  not  one  whose  time  average 
is  equivalent  to  the  mean  realization. 

In  the  analysis  of  nonstationarv  phenomena,  there  are  a  certain  properties  of  the 
resulting  spectrum  which  must  be  identical  to  the  stationary  analog.  These  properties 
include  an  all-positive  spectrum  and  zero  energy  in  the  spectrum  when  the  signal  is  not 
present.  Further  constraints  must  be  applied  to  the  spectral  behavior  along  the  time 
dimension,  a  problem  unique  to  the  anaysis  of  nonstationarv  phenomena.  There  have 
been  many  attempts  to  adequately  model  the  time-varying  behavior  of  nonstationarv 
spectra.  In  general  it  appears  that  each  technique  has  some  advantages  and  disadvan¬ 
tages.  Some  appear  better-suited  to  the  analysis  of  a  certain  class  of  signals.  There  has 
yet  to  be  found  a  completely  satisfactory  description  of  the  time  dependent  spectral  es¬ 
timation  problem. 
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II.  SPECTRAL  ANALYSIS  OF  SIGNALS  WITH  STATIONARY 

CHARACTERISTICS 

The  Wiener- Khinchin  theorem  states  that  if  je(r)  is  a  band  limited,  wide  sense  sta¬ 
tionary  process,  then  the  power  spectral  density  (PSD),  is  related  to  the  autocorrelation 
function  (ACF),  through  a  Fourier  transform 

pjf)  =  P  R^ty***,  (i) 

J-oo 

where  Rxx{r)  is  the  ACF.  Using  the  Fourier  inversion  formula, 

=  [“  (2) 

~oo 

and  evaluating  the  correlation  function  at  lag  zero  results  in  the  average  power  in  the 
process. 

Urn  -jf  J  r  £[  WO  I  !]<*  - 1“  fJJW.  (3) 

where  E  denotes  statistical  averaging  and  2T  represents  the  duration  of  observation. 
Equations  (1)  and  (2)  describe  the  relationship  between  the  time  dependent  data  and  the 
frequency  dependent  power  spectrum.  Equation  (3)  describes  the  relationship  between 
a  signal's  temporal  density  and  its  spectral  density.  The  integrand  to  the  left  in  (3)  re¬ 
presents  the  instantaneous  power  of  the  process.  The  integrand  to  the  right  represents 
the  power  as  a  function  of  frequency.  Both  can  be  considered  power  density  functions 
and  both  are  non-negative  everywhere.  Because  the  ACF  is  always  conjugate  symmetric, 
the  power  spectral  density  is  always  real  and  the  average  power  in  the  process  is  always 
real.  [Ref.  1,2] 
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Equation  (1)  can  be  rewritten, 


PA =  £j 


r-oo 


1 

2r 


x(t)e~jW‘dt 


(4) 


which  implies  knowledge  of  the  signal  for  all  time.  It  is  clear  that  in  order  to  compute 
the  true  ACF  or  PSD,  an  infinite  data  set  is  required.  Knowledge  of  an  infinite  number 
of  realizations  is  also  implied.  These  two  theoretical  constraints  are  impractical.  Real¬ 
istically  one  must  deduce  the  power  spectrum  from  a  single,  finite  duration  realization. 
This  chapter  examines  a  variety  of  estimation  techniques,  noting  the  specific  advantages 
and  limitations  inherent  in  each. 

A.  CLASSICAL  TECHNIQUES 

The  computation  of  a  power  spectrum  from  one,  finite  set  of  data  serves  as  an  esti¬ 
mate  of  the  true  PSD.  One  common  estimation  technique  is  the  periodogram.  It  can 
be  computed  from  the  data  as 


f  x(t)e~i2'Pldt 


Jo 


2 


(5) 


an  estimate  similar  in  form  to  (d),  except  that  limiting  and  statistical  averaging  oper¬ 
ations  have  been  ignored.  This  estimate  has  the  advantage  of  being  both  real  and  posi¬ 
tive.  Examination  of  the  mean  and  variance  of  the  periodogram  spectral  estimate  best 
describes  its  deviation  from  the  true  PSD, 


(*oo 


e[pfM]  = 


sinc[T(f -  a)] Pxx(a )da, 


(6) 


where  it  is  apparent  that  the  periodogram  represents  a  smeared  version  of  the  true  PSD. 
The  smoothing  along  the  frequency  axis  is  caused  by  the  finite  observation  interval. 
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Increasing  the  observation  interval  T  narrows  the  main  lobe  of  the  sine2  function, 
thereby  minimizing  the  smoothing  effect  of  the  convolution  in  (6).  The  variance,  on  the 
other  hand,  is  not  so  accommodating.  [Ref.  3] 

The  variance  of  the  periodogram  spectral  estimate  for  white  Gaussian  noise  is  a 
constant,  with  a  standard  deviation  on  the  order  of  the  mean.  With  so  great  a  variance, 
the  utility  of  the  periodogram  as  described  by  (5)  is  questionable.  One  technique  used 
to  make  this  estimator  more  reliable  is  to  compute  an  average  periodogram.  This  has 
the  effect  of  scaling  down  the  variance  by  the  number  of  terms  in  the  average.  In  prac¬ 
tice,  since  more  than  one  realization  is  rarely  available,  this  amounts  to  segmenting  the 
data  into  shorter  intervals,  which  causes  a  corresponding  increase  in  the  bias  and  loss 
of  resolution  of  the  estimate. 

One  variation  of  the  averaged  periodogram  scheme  requires  a  data  window.  Look¬ 
ing  at  smaller  segments  of  the  data,  the  window  is  applied  in  an  overlapped  fashion  and 
subsequent  processing  results  in  a  series  of  periodograms  that  tend  to  be  correlated  and 
hence  statistically  dependent.  Consequently,  the  actual  reduction  in  variance  will  gen¬ 
erally  be  less  than  the  number  of  terms  in  the  average.  Still  another  variation  requires 
the  data  be  prewhitened.  This  has  the  effect  of  reducing  spectral  bias,  which  is  a  prob¬ 
lem  compounded  by  the  averaging  process.  [Ref.  3] 

Increasing  the  data  length  increases  the  resolution  of  the  periodogram.  An  im¬ 
provement  in  variance  can  be  realized  if,  instead  of  a  periodogram,  one  uses  a 
Blackman-Tukey  spectral  estimator.  This  estimator  is  derived  from  a  biased  ACF  esti¬ 
mate, 


RxxW  =  Tf  J  +  r)di,  (7) 

where  T  is  the  duration  of  the  data  interval.  Taking  the  Fourier  transform  of  a  win¬ 
dowed  version  of  (7), 


hiW= 

(8) 

W[f-  a)Pper(a)da, 

leads  to  an  improvement  in  variance.  The  type  of  window  indicated  by  (8)  should  be 
capable  of  enhancing  the  spectral  characteristics  of  interest  [Ref.  4].  An  unavoidable 
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loss  in  fine  spectral  detail  is  the  price  paid  to  improve  the  variance  of  the  spectral  esti¬ 
mate.  Here  2L  and  T  are  the  durations  of  the  lag  window  and  the  total  data  observation 
respectively.  The  variance  of  the  Blackman-Tukcy  estimator  is 

versus 

(W) 


for  the  periodogram.  [Ref.  3] 

Classical  techniques  are  reliable  but  have  limited  resolution  and/or  a  poor  variance. 
One  assumes  that  the  data  is  zero  before  and  beyond  the  observation  interval.  This  as¬ 
sumption  causes  the  resulting  spectral  estimates  to  deviate  from  the  true  PSD.  Using 
modern  techniques,  this  type  of  error  can  be  minimized.  There  are  however,  other  limi¬ 
tations  to  consider. 

B.  MODERN  TECHNIQUES 

Modern  spectral  estimation  techniques  rely  on  linear  filter  theory.  Parametric 
modeling  provides  the  foundation  for  many  modern  spectral  estimation  procedures. 
Rather  than  assuming  the  data  to  be  zero  beyond  some  arbitrary  interval,  parametric 
models  assume  statistical  knowledge  of  the  underlying  process.  The  data  is  assumed  to 
be  composed  of  sinusoids  in  white  noise.  The  model  should  accurately  estimate  filter 
coefficients  for  some  linear  filter  whose  output,  when  driven  by  white  noise,  is  the 
available  time  data  and  hence  has  a  spectrum  exactly  like  that  of  the  signal  in  question. 
The  power  spectral  density  for  such  an  output  is 

Pxx(J)=\H(J)\Vn.  (11) 

It  is  assumed  that  the  transfer  function,  H(f),  can  be  written  as  a  rational  polynomial 
and  that  the  coefficients  of  the  polynomial  are  the  necessary  filter  coefficients  with  o*  the 
variance  of  the  driving  noise.  Often  the  resonances  of  a  particular  random  process  are 
of  interest.  In  this  case  an  autoregressive  (AR)  model  can  be  sufficient. 

AR  modeling  of  spectra  is  the  most  popular  modeling  method.  By  solving  a  set  of 
linear  equations,  accurate  estimates  of  the  parameters  can  be  determined.  Selection  of 
an  AR  model  results  in  an  all-pole  filter,  hence  the  apparent  high  resolution.  Selection 
of  model  order,  p,  determines  the  number  of  poles  in  the  filter,  which  determines  the 
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number  of  peaks  in  the  spectral  estimate.  Once  the  coefficients  {«(<•)}(„,  are  known,  the 
AR  spectral  estimator  becomes 


ParV)  = 


(12) 


The  statistics  of  the  AR  spectral  estimator  are  difficult  to  determine  in  closed  form. 
Insight  can  be  gained  into  its  reliability  by  considering  the  following  fundamental  as¬ 
sumptions: 

•  the  process  to  be  modeled  is  truly  an  AR  process, 

p 

•  the  process  contains  —  real  sinusoids  or  p  complex  sinusoids,  and 

•  the  coefficients  {«(/')}?,,  are  accurate. 

The  choice  of  model  order  obviously  requires  some  knowledge  of  the  signal's  spectral 
content.  The  number  of  peaks  in  an  AR  spectrum,  for  complex-valued  data,  is  equal  to 
the  model  order,  p.  The  spectrum,  if  the  number  of  sinusoids  truly  present  differs  from 
p  too  much,  can  be  unreliable.  Furthermore,  if  the  signal-to-noise  ratio  is  poor,  the  es¬ 
timate  can  be  of  poor  quality.  [Ref.  3] 

Other  parametric  models  may  be  more  appropriate  to  the  process  under  investi¬ 
gation.  Moving  average  (MA)  models  are  used  when  the  valleys  or  zeros  of  a  spectrum 
are  important.  Unfortunately,  MA  models  require  the  solution  to  a  set  of  nonlinear 
equations.  Although  not  impossible,  it  is  much  more  difficult  to  arrive  at  accurate  pa¬ 
rameter  estimates.  Still  another  alternative  is  to  use  a  combination  of  AR  and  MA 
modeling,  resulting  in  what  is  called  an  ARMA  model. 

The  minimum  variance  spectral  estimator  (Capon's  method)  provides  reasonable 
estimates  without  making  any  assumption  about  the  data  composition  other  than  it 
originates  from  a  wide-sense  stationary  random  process.  The  classical  periodogram 
spectral  estimate  can  be  interpreted  as  a  bank  of  narrowband  filters,  all  with  identical 
passband  characteristics.  The  minimum  variance  method  is  based  on  an  adaptive 
passband  characteristic  with  the  advantage  being  the  ability  to  adjust  sidelobe  levels  and 
hence  minimize  spectral  leakage.  The  resolution  of  this  method  lies  somewhere  between 
that  of  the  periodogram  and  that  enjoyed  by  AR  estimates.  [Ref.  3] 
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III.  SPECTRAL  ANALYSIS  OF  SIGNALS  WITH  DYNAMIC 

CHARACTERISTICS 


A.  ADAPTING  ESTABLISHED  STATIONARY  TECHNIQUES 

*  In  Chapter  II,  classical  spectral  estimators  were  found  to  be  limited  by  the  interval 
of  observation,  limiting  the  spectral  resolution.  A  finite  duration  signal  permits  only  a 
finite  duration  correlation  estimate,  with  increasingly  poor  estimates  away  from  the  zero 
lag.  Modern  spectral  methods  have  apparent  higher  resolution  but  are  sensitive  to  the 
signal-to-noise  ratio  (SNR).  All  the  estimators  discussed  thus  far  assume  the  data  to 
be  at  least  wide-sense  stationary.  The  analysis  of  nonstationary  phenomena  complicates 
the  estimation  problem. 

A  wide-sense  nonstationary  process  is  one  whose  statistics  or  parameters  vary  with 
time  [Ref.  5].  The  actual  nonstationarities  of  a  process  may  include  fluctuating  power 
magnitude,  changing  frequency  content  or  a  combination  thereof.  There  remains  the 
problem  of  a  suitable  representation,  one  which  appropriately  displays  the  fluctuation 
in  time  and  frequency  of  the  instantaneous  energy  in  the  signal.  To  adapt  stationary 
estimation  techniques  to  the  nonstationary  case,  the  concept  of  local  stationary  is  in- 

*  troduced.  A  process  is  considered  locally  stationary  if,  over  a  given  interval  of  time,  the 
process  appears  to  be  stationary.  Determination  of  the  optimum  interval  depends  upon 
the  most  rapidly  fluctuating  nonstationarity  present  in  the  process.  If  the  process  is 
observed  overly  long,  the  temporal  fluctuations  will  be  smeared.  If  the  interval  is  un¬ 
necessarily  short,  spectral  detail  will  be  lost.  A  discussion  of  the  more  popular  methods 
of  nonstationary  estimation  and  their  limitations  follows. 

1.  Short-time  Fourier  Transform 

Short-time  Fourier  analysis  is  a  method  whereby  the  observed  signal  is  seg¬ 
mented  into  a  number  of  shorter  intervals.  Each  segment  is  Fourier  transformed  and 
then  magnitude  squared.  The  resulting  spectra  are  interpreted  as  cross  sections  of  the 
true  instantaneous  spectrum.  These  cross  sections  are  pieced  together  sequentially  to 
form  an  estimate  of  the  true  time-varying  power  spectrum  [Ref.  5,  6].  This  type  of 
spectrum  is  called  a  waterfall  display  by  the  signal  processing  community.  Processing 

* 

sequential  sections  of  data  results  in  a  very  crude  estimate  of  the  nonstationarities  as  a 
function  of  time.  Furthermore,  using  longer  segments  increases  spectral  detail  but  tends 
to  average  or  broaden  the  time-dependent  fluctuations.  An  overly  short  interval  will 
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tend  to  enhance  the  detection  of  the  time  transient  behavior  of  the  process  at  the  ex¬ 
pense  of  spectral  resolution. 

2.  Spectrogram 

A  modified  version  of  the  short-time  Fourier  transform  estimation  technique  is 
the  spectrogram  It  is  related  in  that  it  uses  a  real,  finite  duration,  sliding  window 
(Ref.  7]  centered  at  time  t.  The  spectrogram 


A 

Pspearog. 


(f'l) 


J  .x(t)w(t  ~i)e~^rdx 


(13)  . 


computes  a  classical  estimate  of  the  spectrum  for  each  point  in  time  rather  than  for 
contiguous  blocks  of  data  (see  for  example,  Figure  1).  Similar  to  the  periodogram,  this 
spectral  estimate  is  real  and  positive  everywhere.  The  reliability  of  the  spectrogram 
hinges  on  its  ability  to  represent  the  signal's  energy  in  the  time-frequency  plane.  Refer¬ 
ring  back  to  (3),  any  suitable  time-frequency  representation  should  reduce  to  |jc(/)|j 
when  the  frequency  dependence  is  removed.  Likewise,  the  representation  should  reduce 
to  |A’(/)|J  when  the  time  dependency  is  removed.  Taking  a  look  at  the  spectrogram's 
behavior, 


(a)  amplitude  (b)  contour 

Figure  1.  Behavior  of  the  Spectrogram  using  an  analytic  sinusoid 
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we  note  that  the  estimate  is  smeared  not  only  in  frequency  but  in  time  axis  as  well,  a 
result  of  the  sliding  window  operation.  Another  indicator  of  the  spectrogram's  reliability 
can  be  found  by  removing  both  the  time  and  the  frequency  dependency.  The  equality 
expressed  in  (3)  for  stationary  phenomena  should  have  a  nonstationary  counterpart. 
Therefore  the  volume  under  the  spectrogram  should  equal  the  average  energy  in  the 
signal: 


'oo 

('oo 

*oo 

'oo 

P Specirog.if = 

|  x(t)  I  2vv2(t  —t)chdt.  ( 1 6) 

*'-oc* 

-oo  * 

-oo^ 

— OO 

The  total  signal  energy  requirement  is  satisfied  only  if  the  average  energy  of  the  window 
is  equal  to  unity.  The  spectrogram  represents  a  closer  approximation  of  the  instanta¬ 
neous  energy  changes  relative  to  the  short-time  Fourier  transform;  however,  temporal 
smearing  is  now  unavoidable.  [Ref.  7,  8] 

3.  Correlator/Matched  Filter 

The  spectrogram  estimator  can  be  implemented  as  a  bank  of  bandpass  filters, 


A  „ 

Pcorrifo-1)  =  J  x{i)hfc{t  +  r)dr 


(17) 
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centered  at  each  frequency  of  interest,  f  and  weighted  according  to  some  general  lowpass 
characteristic,  w(  -/)  .  The  output  of  each  bandpass  filter  is  then  magnitude  squared, 
creating  an  instantaneous  energy  representation  of  the  signal  as  a  function  of  time  and 
center  frequency.  Since  the  signal  is  known,  the  general  passband  characteristic  can  be 
replaced  by  a  matched  filter.  Commonly  used  in  a  radar  environment,  each  filter  tests 
for  round  trip  time  delay  and  Doppler  shift.  The  filters  possess  finite  bandwidth  and 
hence  spectral  resolution  finer  than  dictated  by  the  filters  is  impossible.  Similarly,  re¬ 
solution  along  the  time  axis  is  limited  by  the  impulse  response  of  the  filters  [Ref.  5). 

4.  Autoregressive  Modeling 

Autoregressive  modeling  can  be  adapted  to  fit  slowly-varying  spectral  charac¬ 
teristics.  How  slowly  the  frequency  fluctuations  must  occur  depends  on  the  actual 
process  in  question.  In  general,  as  the  signal  rotates  away  from  the  pole  of  the  AR  filter, 
broadening  of  the  peak  results.  As  the  correlation  function  is  time-dependent,  suitable 
accuracy  in  the  AR  coefficients  requires  an  embedded  time  dependence.  This  particular 
modification  is  discussed  in  a  later  section. 

B.  TIME-DEPENDENT  SPECTRAL  TECHNIQUES 

Adapting  stationary  techniques  to  the  nonstationary  case  is  only  marginally  suc¬ 
cessful.  There  have  been  many  attempts  to  describe  the  variation  of  signal  energy,  a 
function  of  both  time  and  frequency,  as  a  multivariable  density.  More  appropriately,  the 
simultaneous  distribution  of  signal  energy  in  time  and  frequency  requires  definition. 
Each  technique  has  its  merits,  and  each  its  own  peculiarities.  In  1966,  a  generalized 
phase-space  distribution  function  was  proposed  [Ref.  9],  which  can  be  used  to  derive 
many  of  the  more  popular  time-frequency  representations.  This  generalized  distribution 
function  is 


moo 

<I>(u,  r).v(q  +  f )/(/,  -  f  )eJ2^!>  -Vt-^dv  dt ,  dr. 


(19) 


The  advantage  in  using  the  generalized  distribution  (19)  lies  in  the  ability  to  define 
properties  belonging  to  all  representations  derived  in  this  manner.  The  distribution  de¬ 
pends  on  the  choice  of  <I)(o,  r),  referred  to  as  the  kernel  of  equation  (19).  An  excellent 
presentation  of  the  relationship  between  particular  time-frequency  distributions  and  (19) 
can  be  found  in  [Ref.  10)  and  is  summarized  below. 
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1.  The  Running  Spectrum 

a.  Using  Only  the  Past  and  Current  Data 

In  1952  C.H.  Page  [Ref.  11]  derived  a  time-frequency  presentation  arguing 
that  one  has  knowledge  of  the  signal  up  to  and  ind’  ding  time  t,  but  its  future  values  arc 
unknown.  He  defined  a  running  transform,  loohing  backwards  over  all  previous  data 
as 


A7 (/)=['  xitJe-M'dti,  (20) 

“00 

where  the  superscript  (-)  indicates  that  the  signal  has  been  observed  over  the  interval 
( -  oo,  /).  By  differentiating  the  squared  magnitude  of  (20)  with  respect  to  time 

?“(/;/)  =  ^- I  A7(/)  |2,  (21) 

the  running  spectrum  suggested  by  Page  results.  Substituting  (20)  into  (21)  and  com¬ 
puting  the  partial  derivative  leads  to  the  following  alternate  form 

=  (22) 

M.H.  Ackroyd  has  argued  [Ref  5]  that  (22)  is  a  finite  duration  approximation  relative 
to  the  physical  measurement  of  a  true  time-varying  energy  distribution.  He  suggests  that 
the  true  time-varying  spectrum  for  a  real-valued  signal  can  be  expressed  as 

*(/)fo[A'(/)e/2c/']  =  x(t)  f  °°  jc(t)  cos [2nj[t  -  z))dr,  (23) 

^ — OO 

the  product  of  the  response  of  a  linear  filter  driven  by  the  signal  and  the  signal  itself. 
The  implementation  suggested  by  (23)  is  shown  in  Figure  2  (a).  It  requires  an  infinitely 
narrow  filter  with  a  noncausal  impulse  response.  To  deal  with  the  causality  issue,  the 
impulse  response  is  modified  by  a  unit  step  function, 


x(t) 


’oo  r  i 

x(t)u(i  -  t)  cos{2nJ{t  -  z))dx  -  x(t)  x(t)  cos(2 nj[t  -  r ))dx 

— OO  OO 


=  x(i)Rel A7(/y2?r/'] 

= y  p~m 


(24) 
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(c)  IPS 

Figure  2.  Time-Frequency  Distribution  Models 
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which  is  a  scaled  version  of  (22).  Figure  2  (b)  suggests  a  practical  means  of  measuring 
a  time-dependent  spectrum.  For  ease  of  comparison.  Figure  2  (c)  shows  the  imple¬ 
mentation  of  a  distribution  discussed  on  pages  14-15.  The  distribution  of  a  signal's  en¬ 
ergy  as  proposed  by  Page  is  limited  in  spectral  resolution;  i.e.,  limited  by  the  finite 
bandwidth  of  the  filter.  It  represents  an  improvement  over  the  spectrogram  which  was 
found  to  be  smeared  in  both  the  time  and  frequency  directions. 

Page's  distribution  can  be  generated  from  the  generalized  equation  (19)  us¬ 
ing  the  kernel  function 

(p(a,T)  =  ^uit|)  (25) 

which  allows  one  to  write  (21)  in  yet  another  form. 
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x  (t)x(l  +  r)e~P"J‘dr  + 
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*0 


=  2Rc[x 


Equation  (26)  can  be  interpreted  as  the  Fourier  transform  of  an  estimate  of  the  true 
time-varying  ACF,  where 

R{[,  t)  =  /(t).v(r  +  t)  for  -  oo  <  r  <  0  ^ 

=  x(i)x  (t  —  t)  for  0  <  r  <  oo. 

The  behavior  of  this  distribution  can  be  seen  in  Figure  3  (a).  As  time  increases,  the 
amplitude  of  the  signal  continually  increases  and  the  true  frequency  location  becomes 
more  localized. 
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0  mxcxi  p$ii 

(a)  amplitude  (b)  contour 

Figure  3.  Behavior  of  Page's  distributions  using  an  analytic  sinusoid 

b.  Using  Only  the  Future  and  Current  Data 

In  1967,  M.J.  Levin  [Ref.  1 2]  extended  the  concept  of  a  running  transform 

to  include 

(28) 

Jt 

where  the  superscript  ( + )  indicates  only  those  data  values  occurring  at  or  later  than  time 
t  are  to  be  considered.  The  equivalent  to  (22)  is 

rVA  —  iW  (29) 

-2Rc[/m7U)elu/ll 

and  the  corresponding  kernel  function  is 

<!>(o,  t)  =  e-'/nu|T|.  (30) 

The  behavior  of  this  future  term  can  be  seen  in  Figure  A.  Maximum  frequency  local¬ 
ization  occurs  early  in  the  distribution,  decreasing  continually  as  time  progresses. 

c.  Using  All  of  the  Data 

Assigning  equal  weight  to  the  past  and  future  terms  Levin  defined  the  in¬ 
stantaneous  power  spectrum  (IPS)  as  the  average  of  the  two  running  spectra, 


ncauDCT  -  mis 
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(a)  amplitude 

Figure  4.  Behavior  of  Levin's  distribution  for  an  analytic  sinusoid 


//>s(/;/)=4[a?(/)  +  .17  (/)] 

=  Re[x\i)e,2nfXPt{f)  +  P7m  1  (3I) 


Figure  2  (c)  suggests  a  method  whereby  the  IPS  may  be  generated.  IPS,  like  Page's 
distribution,  requires  a  noncausal  infinitely  narrow  bandwidth  filter.  Modifying  the  im¬ 
pulse  response  in  order  to  create  a  realizable  filter  (Ref.  13,  14  ]  causes  both  temporal 
and  spectral  smoothing,  [Ref.  15,  pp,  26-28].  The  two  terms  of  IPS  can  be  interpreted 
as  follows/  The  past  term  contains  information  of  the  energy  and  energy  flow  to  create 
the  signal  up  to  time  t.  The  future  term  contains  the  information  about  the  energy  and 
energy  flow  of  the  signal  after  time  t  (Ref.  6]. 

2.  Instantaneous  Power  Spectrum  (IPS) 

IPS  can  be  derived  from  the  generalized  time-frequency  distribution  using  the 
kernel  function: 


<I>(u,  t)  =  cos(7u>t).  (32) 

This  kernel,  cos  rror,  can  be  formed  by  taking  one  half  the  sum  of  the  kernels  for  the  past 
and  future  running  spectra.  Substituting  (32)  into  (19)  and  simplifying  gives 
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(a)  amplitude 

Figure  5.  Behavior  of  IPS  for  an  analytic  sinusoid 


(b)  contour 
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(33) 


It  is  important  to  note  that  the  terms  m  the  integrand  of  (33)  do  not  bear  a  direct  re¬ 
lationship  to  the  past  and  future  spectra  defined  previously  by  (31).  In  this  form,  each 
term  in  the  sum  spans  all  the  data  and  therefore  contains  contributions  from  both  run¬ 
ning  spectra  [Ref.  15].  The  ACF  estimate  as  defined  by  IPS  is 

I<iPs{t > T)  -  4"  ( *(0  x*(l  ~  T)  +  x*(0  •*(*  +  r) )  for  -  oo  <  r  <  oo.  (34) 


By  comparing  Figure  3  and  Figure  4,  the  behavior  of  IPS  as  shown  in  Figure  5,  dem¬ 
onstrates  an  improvement  in  end-point  resolution  where  the  main  ridge  is  most  narrow 
at  the  center  of  the  duration. 

3.  Rihaczek  Distribution  (RD) 

Derived  from  physical  considerations  [Ref.  5,161,  the  complex  energy  or 
Rihaczek  distribution  (RD)  is 
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RD(f,t)  =  x{t)X*  (J)e~'2'f' 


(35) 


Since  the  real  part  of  a  complex  function  is  equal  to  the  real  part  of  the  complex  conju¬ 
gate  of  that  function,  it  is  obvious  the  IPS  as  defined  in  (31)  is  equivalent  (33),  the  real 
part  of  (35).  This  relationship  is  depicted  in  i  igure  2  (c).  As  yet,  there  is  no  satisfactory 
interpretation  of  the  imaginary  part  of  the  Rihaczek  distribution  (ImRD),  although  its 
computation  leads  to  an  increase  in  spectral  localization  over  that  of  IPS  for  certain 
signals.  A  closer  look  at  the  behavior  of  IPS  and  RD  can  be  found  in  Chapter  IV. 

There  exists  a  relationship  between  the  spectrogram  estimate  and  the  RD.  Re¬ 
writing  (13) 
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where  RD(f,i )  is  defined  by  ^35).  Equation  (36)  shows  the  spectrogram  to  be  the  2-D 
convolution  of  the  RD  of  the  data  with  the  RD  of  the  window  function.  [Ref.  6,  15] 

The  complex  energy  distribution  can  be  generated  from  (19)  using  the  kernel 

function 


O(o,  t)  = 


(37) 


Substituting  (37)  into  (36)  results  in 
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Equation  (38)  can  also  interpreted  as  the  Fourier  transform  of  an  ACF  estimate,  where 

R{r,  t)  =  x{t)x\i  -  t)  for  —  OO  <  T  <  OO.  (39) 

Because  RD  is  complex,  this  ACF  estimate  cannot  be  an  even  function  of  the  shift  var¬ 
iable,  t.  This  suggests  that  it  is  the  nonstationarities  of  a  process  which  lead  to  an  ACF 
which  is  partially  odd.  The  behavior  of  this  particular  distribution  is  shown  in  Figues  5 
and  6.  Figure  5  shows  the  behavior  of  IPS  for  an  analytic  sinusoid.  IPS  is  equivalent 
to  the  real  part  of  the  RD.  The  imaginary  part,  shown  in  Figure  6  (a)  and  (b),  demon¬ 
strates  an  improved  sensitivity  to  rapid  changes  in  signal  energy  as  a  function  of  fre¬ 
quency.  The  behavior  of  the  RD  is  discussed  in  more  detail  in  Chapter  IV. 

4.  Wigner-Ville  Distribution  (WD) 

Originally  introduced  by  Wigner  in  a  quantum  mechanical  context  [Ref.  17]  and 
extended  to  signal  analysis  applications  by  Ville  [Ref.  18],  the  Wigner-Ville  distribution 
(WD)  is  another  valid  representation  of  a  signal's  energy  as  a  function  of  both  time  and 
frequency.  The  WD  can  be  generated  from  (19)  using 

Q(o,t)=  1  (40) 

as  the  kernel  function.  Substituting  (40)  into  (19)  results  in 
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(a)  amplitude  of  imaginary  part  (b)  contour  of  imaginary  part 


Figure  6.  Behavior  of  IniRD  for  an  analytic  sinusoid 


+  (41) 

»'-oo 

The  WD  can  be  interpreted  as  the  Fourier  transform  of  an  ACF  estimate  defined  by 

R\vd(1>  t)  =  X* (t  -  Y )  x(t  +  y )  for  -  oo  <  T  <  OO.  (42) 

Since  the  WD  is  always  real,  this  ACF  estimate  possesses  even  symmetry  about  the 
point  of  zero  lag.  An  example  of  the  WD  is  shown  in  Figure  7. 

5,  Time-varying  Autoregressive  Models 

An  appropriate  AR  model  for  time- varying  spectra  contains  an  embedded  time 
dependency.  In  so  doing,  the  model  would  be  able  to  track  a  spectral  peak  minimizing 
the  effects  of  broadening.  Rewriting  equation  (12)  to  reflect  this  time  dependence, 

hm-—} — '■ — -•  («) 

1=0 

where  the  coefficients  are  given  by 
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(a)  amplitude 
Figure  7.  .  Behavior  of  WD  for  an  analytic  sinusoid 

K 

*l(t)  =  Yj  aW  M)-  (44) 

k= 0 

The  {%}£.„  represent  the  weights  of  a  sum  of  ortho-normal  basis  functions.  An  appro¬ 
priate  selection  of  ortho-normal  basis  set  ideally  uses  some  a  priori  knowledge  of  the 
spectrum  under  investigation.  The  time- varying  AR  model  order  is  {K  +  1)P  ,  requiring 
the  estimation  of  ( K  +  1)  times  more  coefficients  relative  to  the  stationary  analog. 
Whether  one  can  accurately  model  the  process  with  a  manageable  number  of  coefficients 
will  depend  on  the  particular  process  at  hand. 


recacNci  •  aiis 


(b) contour 
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IV.  COMPARISON  OF  T-F  DISTRIBUTIONS 

A.  THEORETICAL  RELATIONSHIPS 

Equation  (36)  shows  that  the  spectrogram  can  be  represented  as  the  2-D  convo¬ 
lution  of  the  RD  of  the  signal  with  the  RD  of  the  window.  A  similar  expression  can  be 
derived  relating  the  spectrogram  to  the  WD.  Both  RD  and  WD  can  be  derived  from  the 
generalized  distribution  formula  (19),  as  can  many  other  t-f  representations.  It  turns  out 
that  any  particular  distribution,  C(f,t),  whose  kernel  function  satisfies 

<D(t>,T)a>(  — t>,T)=  1,  (45) 

can  be  related  to  the  spectrogram  in  the  following  way  [Ref.  10], 


/•oo  /* oo 


Spectrog 


m = 


c>,  r)  CJf—  o,  T  -l)do  dr, 


(46) 


where  C,  and  CK  are  the  generalized  distribution  functions  of  the  signal  and  window  re¬ 
spectively.  The  spectrogram  itself,  can  be  represented  through  the  generalized  equation 
using  a  kernel  function 


<l>(o.  t)  =  J' 


w(t  +  f)w*(r  +  f)^>' 


di . 


(47) 


The  right  side  of  the  equality  in  equation  (47)  is  the  ambiguity  function  (AF)  of  the 
window,  vv(r)  [Ref.  6].  For  clarity,  the  AF  is  defined  as  [Ref.  8] 


f  OO 

Xjc>.  ?)  =  x  {t+  r)e  J  Mdl 

J-co 

f  OO 

=  -  y )  +  y  )e~j2rv,'dtl 

— oo  ^ 

r~ 

J*  oo 

X{t^^)x\h-^)e>^dh 

— oo 
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So  far,  we  have  seen  that  a  relationship  exists  between  the  spectrogram  and  certain 
other  distribution  functions  (46),  and  that  the  spectrogram  can  be  generated  using  a 
modified  form  of  the  AF  as  its  kernel.  In  fact,  all  the  distributions  discussed  thus  far  can 
be  related  through  the  AF.  Rewriting  (19)  in  a  slightly  different  form, 


Of,t)  =  f  [  f  4>(o,  r)4r,  +  y )/(/,  -  y  )</2*K  °'~fx)do  dt]  dr 

~oc  oo  — oo  “  ** 

mroo  _  ,  _  n  . 

<P(o,  t)  J  *(/,  +  Y  )x  (h  -  y  1  dt\ 

oo 


(49) 


r  oo 

j  *(b  +  y  )**(h  -  y  )^/2::u,1  dh 


where  the  generalized  distribution  is  shown  to  be  the  2-D  convolution  of  the  double 
Fourier  transform  of  the  kernel  with  the  double  Fourier  transform  of  the  complex  con¬ 
jugate  of  the  AF.  Table  1  lists  the  distributions  discussed  thus  far,  along  with  their  re¬ 
spective  kernels. 


Table  1.  DISTRIBUTIONS  AND  CORRESPONDING  KERNEL  FUNCTIONS 


Distributions 

CM 

Kernel  -  3>(c.  t) 

Spectrogram 

|  fx  x(t)  w(t  —t)  e'jlnftdt  | 2 

Page 

i'w 

M 

e)-'-— 

Levin 

-iu'12 

M 

T 

IPS 

Re  0(r)  X’(J)  e~W] 

COS  7TUT 

RD 

x(t)  X'if)  e~M‘ 

g-jnt 

WD 

1 

B.  GENERAL  PROPERTIES 

Comparing  Figure  1  and  Figures  3  through  7,  it  is  apparent  that  the  particular 
kernel  function  has  a  tremendous  influence  on  the  particular  properties  of  the  resulting 
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distribution.  From  a  spectral  point  of  view,  an  appropriate  t-f  spectral  representation, 
hence  appropriate  kernel  function,  should  ensure  certain  properties.  These  properties 
and  corresponding  kernel  restrictions,  if  any,  are  listed  in  Table  2  below. 
[Ref.  10,  15,  19,  20,  21] 

Of  the  many  distributions  discussed,  none  possess  all  the  desired  characteristics.  The 
choice  of  a  particular  t-f  representation  depends  on  the  application  at  hand.  For  prac¬ 
tical  applications,  the  properties  ascribed  to  the  various  distributions  in  Table  2  must 
be  re-evaluated.  When  using  windowed  data,  the  resulting  WD  is  referred  to  as  the 
pseudo- Wigner  distribution  (PWD).  Similarly  modified  versions  of  IPS  and  RD  are  in¬ 
dicated  by  the  subscript  (y),  where  y  =  *(0  w{t)  is  the  product  of  the  data  with  some 
window  function.  Considering  the  linear  operation  of  filtering,  predictable  distortion  of 
the  signal's  true  spectrum  is  encountered.  In  Table  3  a  summary  of  the  effects  of  four 
linear  operations  as  they  relate  to  time-dependent  spectra  is  given.  [Ref.  10,  15,  19] 

C.  RELATIVE  PERFORMANCE 

In  the  previous  section  it  was  shown  that  the  spectrogram  is  related  to  certain  gen¬ 
eralized  distributions  through  a  2-D  convolution.  It  was  further  shown  that  any  gener¬ 
alized  distribution  is  the  2-D  convolution  of  the  double  Fourier  transform  of  the  kernel 
function  with  the  double  Fourier  transform  of  the  complex  conjugate  of  the  ambiguity 
function  (AF).  All  these  inter-relationships  are  interesting.  By  looking  at  the  relative 
performance  of  certain  t-f  distributions,  insight  into  the  benefits  and  disadvantages 
characteristic  of  a  specific  representation  are  easily  seen.  The  distributions  compared  in 
this  section  are: 

•  Spectrogram 

•  IPS, 

•  ImRD ,,  the  imaginary  part  of  RDy 

•  2  linear,  magnitude  combinations  of  the  real  and  imaginarv  parts  of 

RDy 

•  PWD, 

where  the  subscript  (y)  indicates  the  use  of  windowed  data. 

1.  Experimental  Analysis 

Seven  test  signals,  in  the  absence  of  noise,  where  considered.  What  follows  is 
a  brief  description  of  the  true  time  frequency  behavior  of  the  analytic  signals  and  a  list 
of  the  particular  processing  parameters  of  the  digital  implementations.  Unless  otherwise 
specified,  all  data  is  128  points  in  duration,  using  a  1024  point  fFT  algorithm.  Two 
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Table  2.  GENERAL  PROPERTIES  OF  T-F  DISTRIBUTIONS 


Properties 

Kerne!  constraints 

Distributions 

x(l)  -  C(f,t) 

1 

2 

5 

6 

Zero  energy 

x(tB)  =  0  -+  Of,  t0)  =  0 

1 

B 

fl 

1 

X(f0)  =  0  -  Cfo.  t)  =  0 

1 

B 

a 

a 

fl 

1 

Time  shift 

x(t  -  t0)  ->  C(f,t-  f0) 

0)(o.  t)  must  be  independent 
of  absolute  time  and  fre¬ 
quency 

I 

1 

1 

S 

1 

1 

Freq.  shift 

x(t)  eW  -  Cf-f.  t) 

Positive 

Off )  1>  0  Vft 

F,,,[%t)]  2;  0 

B 

fl 

1 

■ 

1 

1 

Real 

C(f,t)  =  C(  —ft) 

fl 

fl 

fl 

fl 

1 

fl 

Marginal  in  t 

/"  C{f,t)df  =  |-x(r)  |2 

d>(o,  0)  =  1  Vd 

1 

s 

a 

fl 

fl 

Marginal  in  f 

/“  C{f,t)di  =  1  A'(/)  | 

*'  -oo 

t)  =  1  Vt 

1 

fl 

X 

a 

fl 

fl 

Group  delay 

/“  t  c(/;out 

-OO  rp 

lA'COP  "  1 

d>(0,  r)  =  0  Vt 
and 

-T|-d>(y,  t)|c=0  =  0 

1 

1 

Instantaneous^ 

frequency 

f-JC(f,t)df 

WOP  ~f‘ 

O(u,0)  =  0  Vi; 
and 

-jf  T)  1  r=0  =  0 

1 

1 

1 

1 

Legend: 

1  =  Spectrogram,  2  =  Page,  3  =  Levin, 

4  =  IPS,  5  =  RD,  ,  6  =  WD 

1  Analytic  signals  only 
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Table  3.  LINEAR  OPERATIONS 


signal 

Fourier 

transform 

Distributions 

Ideal 

X(l) 

w 

WDx{f,t) 

time  shift 

x(t  -  to) 

X(f)e~J2'f'o 

- 10) 

HDjffJ  -  k) 

WD,{f,t-t0) 

modu¬ 

lation 

x  (t)ej2n'o' 

W-fo) 

window¬ 

ing 

x(t)  w(t) 

A 

•  s 

S 

•*< 

Re  lRDx(f,t)'RDwm 

RD,{f,t)*RDw{f,t) 

WDx{f,iyn-Dj/,t ) 

filtering 

x(t)  *  h(t) 

m  m 

Re  tRDtWRDmVj]] 

RDx{f,t)*RDw{f,t) 

l  VDx(f,ty\VD>(f,l) 

different  Hamming  window  functions  are  used  depending  on  the  type  of  spectra  ana¬ 
lyzed,  stationary,  noil-stationary,  or  a  combination  of  both.  Except  for  the  spectrogram, 
all  distributions  have  been  smoothed  in  the  time  direction  using  a  5-cell  box-car  average, 
centered  at  the  time  of  interest.  The  test  signals  used  were: 

1.  Single  component.,  analytic  sinusoid  computed  as  where  1  <  n  <  128,  using 

a  127  point  Hamming  window, 

2.  Two  component.,  analytic  signal  computed  as  eJ1”2M'n  +  ej2'¥x",  where 
1  ^  n  <  128,  using  a  127  point  Hamming  window, 

3.  Single  component,  analytic  linearly  chirped  signal  computed  as 

f!,i (5-0_10<T3r>)  ,  where  1  <  n  <  128,  using  a  55  point  Hamming  window, 

4.  Two  parallel.,  analytic.,  linearly  chirped  signals  computed  as 

e'2,r (i  0  * ,0( TIT ’)  +  ^2'rlfr  (1£  0 "  ,0<Tfr)),  where  !  <  h  <1 128,  using  a  55  point  Hamming 
window, 

5.  Single  component,  analytic,  quadratically  chirped  signal  computed  as 

e'-’'— O'0-11*  ,  where  1  <  /?  <  128,  using  a  55  point  Hamming  window, 

6.  Multi-component  signal  comprised  of  a  stationary  sinusoid,  a  linearly  chirped  and 
a  quadratically  chirped  component  computed  as 

4-  c,'*,T2r(2$-0_r<Tfr))  +  e'J,T»(J-0"WllM>2)  ,  where  1  <  «  <  128,  using  a  55 
point  Hamming  window, 

7.  Frequency  shift  keyed  (FSK)  signal  computed  as  for  l<n<24  anci 

56  <n<  128,  and  for  24  £n<  56,  using  a  55  point  Hamming  window. 
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2.  Highlights  of  the  Analysis 

To  judge  the  accuracy  of  component  placement  along  the  frequency  axis,  the 
performance  of  IPS y  is  compared  to  that  of  the  spectrogram  in  Figure  8.  Using  a  sta¬ 
tionary,  single  component,  analytic  signal  both  spectra  have  good  end-point  resolution. 


C  racauDCl  *»i*  «  .  c  racocxi  a*i* 


(a)  IPS y  (b)  Spectrogram 

Figure  8.  Amplitude  plots  of  a  stationary'  analytic  sinusoid 

Where  the  spectrogram  presents  a  stationary  spectral  ridge  with  increasing  amplitude 
near  the  center  of  the  t-f  plane,  IPSy  shows  a  fairly  constant  amplitude,  stationary 
spectral  ridge  which  is  wider  by  comparison.  Equation  (36)  describes  the  spectrogram 
as  a  2-D  convolution  of  the  RD  of  the  signal  and  the  RD  of  the  window.  The  apparent 
superior  resolution  of  the  spectrogram  over  IPSy,  the  real  part  of  the  RD,  indicates  that 
the  imaginary  part  of  the  RD  contains  some  important  spectral  information.  The  im¬ 
aginary  part  is  formed  as 


1*00 


ImRDif't)  =  f 


( x(t )  x  (/  —  t)  —  ^*(r)  x(r  +  T))e~J1^'dT. 


(50) 


Using  the  kernel  function 


<t>(o,  t)  —j  sin  7tot,  (51) 

equation  (50)  can  be  derived  from  (19).  It  is  obvious  that  the  kernel  for  the  RD  is  the 
sum  of  the  kernels  for  IPS  and  ImRD.  As  to  the  behavior  of  ImRDy,  a  windowed  version 
of  (50),  demonstrates  an  improved  sensitivity  to  spectral  change  relative  to  PWD  and 
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(a)  Amplitude  (b)  Contour 

Figure  9.  IntRDy  for  a  stationary,  analytic  sinusoid 

IPS)  which  is  shown  in  Figure  9.  In  particular,  using  a  single-component,  linearly- 
chirped  signal,  the  zero  crossings  of  lmRDy  occur  near  the  true  frequency  location.  In 
Figure  10,  the  relative  performance  between  PWD,  IPSy  and  lmRDy  for  a  single¬ 
component,  linearly  chirped,  analytic  signal  can  be  compared.  One  can  see  that  ImRDy 
has  an  improvement  in  resolution  over  that  of  IPS,  comparable  to  that  achieved  by 
PWD.  Where  both  IPSy  and  PWD  rely  on  the  sharpness  of  the  spectral  ridge  for  resol¬ 
ution,  IniRDy  relies  on  the  zero  crossings. 

The  lmRDy  has  large  values  near  the  time  of  frequency  change,  with  the  ampli¬ 
tude  approaching  zero  otherwise.  This  behavior  is  clearly  demonstrated  in  Figure  9 
where  a  single  component  stationary  signal  is  shown  to  have  a  nonzero  imaginary 
spectra,  greatest  in  amplitude  near  the  beginning  and  end  of  the  life  of  the  signal.  In 
Figure  8  (a),  IPSy  presents  a  spectral  ridge  which,  although  rounded  at  the  endpoints, 
grows  increasingly  more  narrow  toward  the  center  of  time.  By  forming  a  linear  combi¬ 
nation  of  IPSy  and  lmRDy,  an  overall  improvement  in  spectral  resolution  for  some 
signals  can  be  achieved.  One  such  linear  combination  is 

I  RDy(f,t)  | 2  =  |  IPSy(f,i)  | 2  +  |  ImRDJf,t)  | 2.  (52) 

Another  combination  can  be  formed  taking  the  difference  of  magnitudes, 


\RDym2-=\IPSym7  -  \ImRDJf,t)\\ 


(53) 


nCOUCNCT  •  MIS 

(c)  Contour  plot  of  ImRD f 


Figure  10.  Contour  plots  of  a  single-component  linear  chirp 

Equation  (52)  results  in  a  nonnegative  spectrum.  Equation  (53),  unfortunately,  can  have 
negative  values. 

For  stationary  data,  j  LDf  | 1  represents  an  improvement  in  the  end-point  resol¬ 
ution  relative  to  IPS ,  as  can  be  seen  in  Figure  11.  When  two  stationary  components 
are  present  IPSt,  ImRD \  and  both  linear  combinations  produce  modulated  spectra.  The 
modulation  effect  is  related  to  the  difference  frequency.  PWD  produces  a  spectrum 
having  cross-terms  midway,  between  the  true  components.  The  cross-terms  are  also  re¬ 
lated  to  the  difference  frequency.  The  relative  behavior  of  IPSy,  \RD,\l  and  PWD  for 
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(a)  Contour  plot  of  IPS t  (b)  Contour  plot  of  |  RD>  | J 

Figure  1 1.  Contour  plots  of  a  single-component,  stationary  analytic  sinusoid 

a  two-component  analytic  signal  is  shown  in  Figure  12.  Using  128  point  data  sets  at  a 
sampling  frequency  of  128,  classical  analysis  predicts  component  resolution  where  the 
separation  is  at  least  0.89  Hz  (Ref.  4].  IPS,  is  capable  of  resolving  two  components 
separated  by  at  least  1.25  Hz.  and  PWD  are  able  to  resolve  two  components 

separated  by  at  least  0,6  Ilz,  an  improvement  over  classical  analysis.  Figure  13  is  a 
graph  depicting  the  relative  location  of  the  spectral  peaks  as  a  function  of  frequency  and 
time.  IPSy,  after  an  initial  settling,  consistently  places  the  spectral  peaks  demonstrating 
a  bias  which  is  not  symmetric.  (/?/),  |i,  similar  to  IPSy,  consistently  places  the  spectral 
peaks  but  with  a  symmetric  bias.  PWD  succeeds  in  making  the  nearest  approximation 
to  the  true  component  locations  demonstrating  a  bias  which  is  not  symmetric.  The 
placement  of  these  spectral  peaks  does  not  appear  to  settle  at  one  location  as  appears 
to  be  the  case  for  the  Rihaczek-derived  distributions;  however  in  a  'nean-squared  error 
comparison,  PWD  appears  to  be  superior. 

The  ImRDy  is  also  capable  of  resolving  two  closely-spaced,  narrow-band,  sta¬ 
tionary  components.  Instead  of  searching  for  spectral  peaks,  ImRDy  characteristically 
detects  the  zero  crossings  which  yield  information  in  this  spectrum.  In  Figure  14,  the 
behavior  of  ImRDy  using  two  closely-spaced  stationary  sinusoids  can  be  seen. 

To  study  the  behavior  of  any  spectral  estimator  of  nonstationary  phenomena 
we  begin  by  considering  a  single-component,  linearly  chirped,  analytic  test  signal.  The 
relative  periormance  of  IPSy,  ImRDy  and  PWD  was  discussed  previously,  see  also 
Figure  10.  The  behavior  of  the  two  linear  magnitude  combinations  can  be  seen  in  Fig- 
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(c)  Amplitude  plot  of  |  RDy\i 

Figure  12.  Amplitude  plots  of  a  tuo  component  stationary,  analytic  sinusoid 

urc  15.  Forming  the  difference,  If?/),!!  creates  a  spectral  ridge  comparable  in  width 
to  IPSy  but  with  a  better  defined  peak.  Detection  using  \RDy\i  should  be  thus  be  im¬ 
proved.  Forming  the  sum,  1  RDyl1  completely  resolves  IPS  into  its  past  and  future  terms, 
the  distributions  defined  by  Page  and  Levin.  Component  location  using  1  /?/),! 2  requires 
detection  of  the  minimum  occurring  between  the  two  resolved  ridges,  making  detection 
using  \RDy\2  questionable.  In  the  presence  of  noise,  the  separation  IPS,  into  compo¬ 
nent  parts  will  lead  to  difficulty  in  interpretation.  The  ability  of  lPSy  to  properly  locate 
the  instantaneous  frequency  for  a  linear  chirped  signal  is  shown  in  Figure  16.  The  lo- 
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frequency 

Figure  13.  Graph  depicting  accuracy  of  stationary  component  placement 

cation  of  the  peaks  in  IPS>  coincide  with  the  true  instantaneous  spectral  locations  for  a 
slow  chirp.  Eatly  in  the  spectrum,  when  the  future  term  in  dominant,  IPS ,  tends  to  place 
the  instantaneous  frequency  higher  than  truth.  On  the  other  hand,  late  in  tnc  spectra 
when  the  past  term  is  dominant,  JPSy  tends  to  place  the  instantaneous  frequency  lower 
than  truth.  Doubling  the  chirp  rate  results  in  a  greater  frequency  ambiguity  as  seen  by 
the  two  terms  which  compose  IPS,.  Looking  at  Figure  16,  this  is  demonstrated  by  the 
apparent  random  placement  by  I PSy  of  the  spectra!  peaks,  neither  the  future  nor  past 
term  seems  to  be  favored  as  the  maximum  peak  location. 

Next,  a  test  signal  composed  of  two,  parallel,  linearly  chirped,  analytic  signals 
is  considered.  The  spectra  resulting  from  IPSy  and  PVVD  is  shown  in  Figure  17.  Similar 


Figure  14.  IniRDy  for  a  tuo-coinponent  stationary,  analytic  sinu..  .d 
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(a)  Amplitude  plot  of  (b)  Amplitude  plot  of  \RDy\2 

Figure  IS.  Relative  performance  for  a  single-component  linear  chirp 

to  the  stationary,  tvvo-componcnt  case,  modulation  and  cross-terms  fluctuating  at  the 
difference  frequency  are  present  in  IPSy  and  PWD,  respectively.  Neither  of  the  two  lin¬ 
ear  combinations  of  magnitudes  of  the  component  parts  of  the  Rihaczek  distribution, 
was  able  to  improve  on  the  resolution  achieved  by  IPSy.  Moving  from  the  peaks  to  the 
valleys,  hnRDy  can  be  seen  to  discriminate  between  the  two  chirps  as  a  function  of  the 
zero  crossings.  This  is  demonstrated  in  Figure  18  where  a  pattern  of  zero  crossings  al¬ 
lows  the  eye  to  discern  the  two  components. 


Figure  16.  Graph  depicting  accuracy  of  IPS  in  locating  the  instantaneous  frequency 


32 


C  r*w**ct  mis  •  0  r*tx«»cr  mis 


(a)  Amplitude  plot  of  IPS ,  (b)  Amplitude  plot  of  PWD 

Figure  17.  Behavior  of  IPSy  and  PWD  for  2  parallel  linear  chirps 

Examination  of  the  spectra  for  more  complicated  signal  suggests  a  relative 
ranking  in  terms  of  resolution  among  the  different  power  distributions  discussed  thus  far. 
The  test  signal  is  composed  of  a  high  frequency  stationary  component  and  two  chirped 
components,  one  a  linear  chirp  and  the  other  a  quadratic  chirp.  Considering  first  those 
spectra  which  estimate  frequency  location  as  the  point  of  maximum  power,  PWD 
produces  the  most  narrow  ridges;  see  Figure  19  (a).  It  su(Fers  from  poor  end-point  re¬ 
solution  and  spectral  cross-terms.  \RDy\l,  seen  in  Figure  19  (b),  produces  well-defined 
periodic  peaks  along  the  instantaneous  frequency  path  of  the  nonstationarv  compo¬ 
nents.  The  stationary  component  ridge  is  broadened,  with  the  modulation  effect  most 
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Figure  18.  ImRDy  of  tuo  parallel  linear  chirps 
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(c)  ltnRDy 

Figure  19.  A  combination  of  stationary  and  nonstationary  spectral  components 

apparent  near  the  end-points.  I RD}\2  would  be  diilicult  to  interpret  as  it  resolves  the 
past  and  future  terms  for  the  nonstationary  components.  Looking  for  the  zeros,  hnRDy 
accurately  describes  the  location  of  the  nonstationary  components,  but  provides  little 
information  during  parts  of  the  stationary  signals  existence.  Again,  Figure  19  (c)  re¬ 
quires  pattern  recognition  to  be  able  to  discern  the  individual,  dynamic  components. 

Considering  a  pulsed  spectra,  such  as  in  FSK  modulation,  \RDy\2  presents  a 
narrow  ridge  with  good  resolution  throughout  the  duration  of  each  pulse.  '1  he  ridge 
width  for  each  pulse  in  the  PWD  spectral  description  is  dependent  upon  the  pulse  du- 
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(a)  \RDy\} 

Figure  20.  Contour  plots  for  a  complex  anlaytic 

ration.  The  longest  pulse  showing  the  most  narrow  ridge,  one  which  is  slightly  more 
narrow  than  what  is  found  with  \RD}\}.  These  two  spectra  are  shown  in  Figure  20. 
PWD  has  a  slow  build  up  and  decay  at  the  ends  of  each  pulse.  There  appears  to  be  a 
ttadc-olf  between  the  width  of  a  spectral  tidge  and  end-print  resolution. 

3.  Test  Case  Results 

Figures  of  the  spectral  dist.ibutions  resulting  from  the  six  t-f  representations  are 
contained  in  this  section.  They  arc  oi deied  by  the  type  of  test  signal  under  analysis. 
a.  Single-component,  Analytic  Sinusoid 

Using  the  spectrogram  as  the  reference,  the  resolution  ability  of  five  addi¬ 
tional  t-f  distributions  is  compared  in  Figure  21  -  Figure  24.  PWD  presents  a  well- 
defined  spectral  ridge  near  mid-plane,  but  sulfers  from  a  sluggish  built-up  and  decay.  IPS} 
presents  a  wider  main-lobe  rclathe  to  PWD  which  is  compensated  somewhat  by  an  in¬ 
crease  in  end-point  resolution.  \RDy\2  provides  the  best  end-point  resolution,  main¬ 
taining  constant  amplitude  throughout  the  plane.  |  RDy\i  proves  to  be  more  sluggish 
than  PWD  and  ImRD}  demonstrates  an  improved  response  near  the  time  of  spectral 
change. 
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(b)  PWD 
FSK  signal 
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(a)  Spectrogram 


TREoxKcr  -  axis 


TRCoueNCT  -  axis 

(b)  IPS, 


(c)  PWD 

figure  22.  Test  signal  1:  contour  plots  for  Spectrogram,  IPS,  and  PWD 
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(c)  ImRDy 

Figure  23.  Test  signal  1:  amplitude  plots  for  \RDyl1,  and  ImRDy 
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b.  Two-component ,  Analytic  Signal 

The  spectrogram  and  all  five  t-f  distributions  shown  in  Figure  25  through 
Figure  28  display  a  distorted  spectrum.  This  distortion  exists  for  a  minimal  time  in  the 
spectrogram  estimate.  In  the  case  of  IPSy  and  PWD,  the  distortion  is  related  to  the 
difference  frequency.  IPSy  shows  modulation  of  each  component;  PWD  contains  addi¬ 
tional  peaks  oscillating  at  the  difference  frequency.  The  modulation  effect  characteristic 
of  Rihaczek-derived  distributions  makes  the  spectra  using  1  RDy  \ 2  and  ImRDy  extremely 
difficult  to  interpret.  |  RDy  j  *  suffers  also  from  the  modulation  effect;  however,  this 
magnitude  combination  enhances  resolution  relative  to  IPSy ,  making  it  easier  to  detect 
the  presence  of  two  components. 

c.  Single-component,  Analytic,  Linearly  Chirped  Signal 

The  inadequacy  of  assuming  local  stationarity  is  clearly  demonstrated  in 
Figure  29  (a)  and  Figure  30  (a)  for  the  spectrogram,  where  the  slope  of  the  instantane¬ 
ous  frequency  line  is  distorted  and  broadened  near  the  end-points.  Both  IPSy  and  PWD 
shown  in  Figure  29  (b)  and  (c),  and  Figure  30  (b)  and  (c)  maintain  a  better  approxi¬ 
mation  to  the  instantaneous  frequency  slope  neat  the  end-points.  Characteristically, 
IPSy  presents  a  broadened  spectral  ridge  whereas  PWD  deca>  s  to  zei  o  at  the  start  and 
stop  of  the  chirp.  Looking  at  the  spectrum  created  by  \RDy\2,  Figure  31  (a)  and 
Figure  32  (a),  the  future  and  past  terms  which  make  up  the  RD  are  clearly  visible.,  Al¬ 
though  | RDy\2  is  an  all-positive  distribution  possessing  many  desirable  properties,  this 
characteristic  resolution  quality  makes  it  difficult  to  interpret  more  complicated  spectra. 
The  two  remaining  Rihaczek-related  distributions  represent  improvements  over  the 
ability  of  IPSy  to  pinpoint  the  instantaneous  frequency  as  a  function  of  time.  |  RDy  \  i 
in  Figure  31  (b)  and  in  Figure  32  (b)  shows  a  spectral  ridge  more  narrow  than  that  for 
PWD.  ImRDy,  while  accurately  locating  the  instantaneous  frequency,  requires  the  de¬ 
tection  of  zero  crossings.  For  particular  detection  schemes  information  thus  provided 
may  be  appropriate. 

d.  Two  Parallel,  Analytic,  Linearly  Chirped  Signals 

For  this  test  signal,  the  spectrogram  spectral  estimates  shown  in  Figure  33 
(a)  and  Figure  34  (a)  is  unacceptable.  The  comments  made  previously  concerning  the 
distortion  in  the  spectra  when  when  two  closely  spaced,  parallel,  stationary  components 
can  also  be  applied  to  this  nonstationary  case  (see  Figure  33  -  Figure  36).  In  the  case 
of  PWD,  not  only  are  cross-terms  oscillating  between  the  true  components  present  but 
the  spectral  ridges  themselves  show  the  effects  of  modulation.  \RDy\2  and  ImRDy  do 
not  show  promise  as  estimation  tools  for  closeh -spaced  frequency  components.  It  is 
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interesting  to  note  however,  that  1 | 2  resolve  the  past  and  future  terms  for  each 
component  similar  to  the  behavior  shown  in  Figure  31  (a).  \RDy\l  shows  very  sharp 
peaks  along  the  modulated  instantaneous  frequency  lines  of  the  two  linear  chirps.  Be¬ 
cause  these  peaks  are  the  largest  peaks  in  the  plane,  1  RI)y\l  from  a  practical  view  point, 
appears  to  be  more  suited  toward  the  analysis  of  this  type  oi  signal  than  PWD  uhcie  the 
cross  terms  periodically  represent  the  dominant  peak. 
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(c)  PWD 


Figure  26.  Test  signal  2:  contour  plots  for  Spectrogram,  IPSy  and  PWD 


43 


(c)  ImRD, 

Figure  31.  Test  signal  3:  amplitude  plots  for  1 7? | 2,  |  1 1  and  ImRDy 
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(c)  ImRD, 

Figure  32.  Test  signal  3:  contour  plots  for  \RDy\2,  | RDy\l  anil  IiuRDy 
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(b)  IPS, 


plots  for  Spectrogram,  IPS y  and  PWD 
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(c)  PWD 

figure  34.  Test  signal  4:  contour  plots  for  Spectrogram,  IPSy  and  PWD 
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(c)  IniRD y 

Figure  36.  Test  signal  4:  contour  plots  for  1  RDy\\  I RDy\i  and  hnRDy 
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e.  Single-component,  Analytic,  Qnadratically  Chirped  Signal 

The  ability  of  the  spectrogram  and  these  5  t-f  distributions  to  accurately 
display  more  rapid  spectral  dynamics  can  be  compared  in  Figure  37  to  Figure  40.  Not 
surprisingly,  the  spectrogram  presents  a  poor  estimate.  IPSy  demonstrates  a  spectral 
ridge  which  tends  to  broaden  as  a  function  of  time,  making  it  difficult  to  ascertain  the 
actual  instantaneous  frequency  curve.  PWD,  although  zero  at  the  end-points,  tracks  the 
chirp  closely  presenting  a  ridge  along  the  line  of  instantaneous  frequency  as  narrow  as 
that  found  for  the  linear  chirp.  |/?Z)j.|l  appears  to  be  provide  the  most  narrow  ridge. 
Both  PWD  and  show  amplitude  modulation  along  the  peak  of  the  curve.  As  to 

the  remaining  spectra,  1 RDJ2  and  lmRDy  behave  in  a  manner  similar  to  the  case  of  the 
linear  chirp  and  thus  do  m  appear  particularly  suited  to  this  class  of  signal. 

f.  Multi-component  Analytic  Signal 

How  the  various  spectral  estimation  techniques  perform  when  confronted 
with  a  mixture  of  stationary  and  nonstationary  dynamics  is  demonstrated  in  Figure  41 
through  Figure  44.  This  test  signal  is  composed  a  high  frequency  stationary  component 
and  two  chirped  components,  a  linear  chirp  and  a  quadratic  chirp.  Their  relative  per¬ 
formance  suggests  a  ranking  in  terms  of  desirability  as  an  estimator  of  continuously 
changing  spectral  components.  The  spectrogram  is  predictably,  the  worst,  followed  by 
|  /?Z>V  | 2  and  IniRDy  in  ascending  order  of  desirability.  Only  IPSy  ,  and  PWD 

present  spectra  which  closely  resemble  the  true  signal  components.  All  three  show  a 
broadened  ridge  for  the  stationary  component  relative  to  the  single  stationary  compo¬ 
nent  case  examined  previously.  This  results  from  using  a  shorter  window.  As  expected, 
| RDy\l  sharpens  the  modulation  peaks  found  in  its  own  spectrum  and  that  of  IPSy.  A 
comparison  of  and  PWD  is  also  quite  characteristic.  The  price  paid  for  in¬ 

creased  end-point  resolution  and  elimination  of  the  spectral  cross-terms  is  a  broadening 
of  the  spectral  ridge  and  the  appearance  of  modulation  along  the  crest  of  the  spectral 
ridges. 

g.  Complex  FSK  signal 

In  contrast  to  the  continuously-varying  frequency  dynamics  of  the  previous 
test  case,  a  very  different  order  of  desirability  is  suggested  in  the  case  of  pulsed  spectral 
dynamics.  Comparison  of  Figure  45  -  Figure  48  shows  |  RDy  \ 2  to  possess  superior 
end-point  resolution  ability  coupled  with  a  rapid  build-up  and  decay  of  the  spectral  ridge 
reiathe  to  the  true  pulse  dynamics.  IPSy  and  j  RDy\ i  produce  similai  jvpectra;  however, 
the  end-point  build  up  is  slower.  JmRDy,  with  its  heightened  sensitivity  to  detect  spectral 
change  appears  to  be  an  excellent  indicator  of  the  time  and  location  of  frequency 
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(c)  PWD 

Figure  37.  Test  signal  5:  amplitude  plots  for  Spectrogram,  IPSy  and  PWD 


change.  PWD  presents  a  sluggish  build-up  and  decay,  with  the  width  of  each  spectral 
ridge  dependent  on  the  duration  of  the  pulse.  Furthermore,  cross-terms  can  be  seen 
between  the  two  higher  frequency  pulses. 
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(c)  PWD 

Figure  38.  Test  signal  5:  contour  plots  for  Spectrogram,  lPSy  and  PV\  D 
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(c)  PWD 

Figure  41.  Test  signal  6:  amplitude  plots  for  Spectrogram,  IPSy  and  PWD 
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figure  42.  Test  signal  6:  contour  plots  for  Spectrogram,  lPSy  and  PWD 
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Figure  44.  Test  signal  6:  contour  plots  for  |  RD ( |:,  |  RDy  |i  and  ImRDy 
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(c)  PWD 

Figure  45.  Test  signal  7:  amplitude  plots  for  Spectrogram,  IPSy  and  PWD 
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(c)  ImRD, 

Figure  47.  Test  signal  7:  amplitude  plots  for  | RDy\\  | /?£>,, 1 1  and  ImRDy 
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(c)  lmRDy 

Figure  48.  Test  signal  7:  contour  plots  for  |  I  RDr  |L  and  ImRDy 
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D.  ALTERNATE  METHODS  OF  COMPUTING  IPS 

Equation  (33),  the  defining  equation  for  IPS,  can  be  rewritten  as  the  Fourier  trans¬ 
form  of  an  ACF  estimate  (34).  namely 


/*oo 


IPS(f,t) 


RiPs(r,l)e~J2'f‘  dr 


(54) 


(xf;)x  (t  —  t)  +  x  {i)x(t  +  r))e  dr. 


IPS  for  finite-duration  discrete  signals  is 


IPS{6.n)-—-‘  (x(n)x  (n  —  k)  +  x’(n)x{n  +  k))e~Jdk,  (55) 

fc=— oo 

where  the  signal  sequence  x{n)  is  finite  and  zero  outside  the  known  samples,  and  AT  is 
a  constant.  Equation  (55)  can  be  expressed  as  [Ref.  15] 

lPS(d,,t)  -  |.v(«)l2=  U’(0)|2  -  \X2(6) I2,  (56) 

where 


oo 

m  =  X  x{r)e-JOr 

oc 

xm=  X  x^e~j0r- 

r=— oo 
r*n 


(57) 


For  use  later  in  the  derivation,  let  the  Fourier  transform  of  the  point  of  interest  be 

IW)|J-W»)|!. 

By  neglecting  |  *(.-/)  | 2  in  (56)  a  modified  version  of  IPS  is  defined.  The  behavior  of  this 
modified  IPS  for  a  single-component,  anal} tic  sinusoid  is  shown  in  Figure  49  (a). 
Comparing  this  with  an  unwindowed  version  of  IPS  for  the  identical  data  sequence  in 
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Figure  49.  IPS 


(a)  Modified  IPS 


(b)  IPS 


Figure  49  (b),  shows  the  two  methods  of  computing  IPS  to  be  very  similar.  The  greatest 

I 

difference  can  be  seen  when  comparing  maximum  amplitudes.  Modified  IPS  is  on  the 
order  of  10-5  and  IPS  is  on  the  order  of  10.  Both  methods  were  computed  using  a  five- 
cell-box-car  averaging  procedure  along  the  t  axis. 

The  question  arises  concerning  the  implementation  of  modified  IPS  using  a  window 
function.  One  method  of  windowing  can  be  written  as 


PC 

IPS(9,n)  =  -~-  (.v(/j)jc*(h  -  A)  +  x* (n)x(n  +  k))w(k)  e~J9k .  (59) 


k=— oo 


Applying  the  definitions  in  (57), 


7PS(0,H)  =  ~2-{A:(«)e";(?''(/(0)*IF(0))  +  x'(n)e,0tt(X(6)*  IF(0))J 


AT 


AT 

2 


{D(d){x’(6)*}V(d))  +  D\d){X{9)HVm} 

oo 

jfl(0)  ^  X'(p)W(B  - p )  +  D\0)  J  X{p)lV{0  -p)j, 


p=-oo 


(60) 


where  *  denotes  convolution  in  the  frequency  variable.  Moving  everything  inside  the 
summation  sign  and  substituting  X[{p)  +  D‘(p)  for  X’{p)  and  A'*(0)  -  A"(0)  for  D'(0 )  gives 


oo 


IPS(d,n)  =  ^f  )  (A>)  +  Z)*(^))D(0)IF(0-P) 


Z-J 

p=— oo 


V  (a-* («)  -  XwM'w»n«-rt 


P=— oo 


-¥•  )  [W)4w  +  *W>’«]w,(«-f) 


OO 


(61) 


Z_J 

p=— oo 


+ 


Ar 


oo 

^  [A'*(0)*(p)  -  Aj(0)A'O»)]H'(0-/»). 


/?=- OO 


A  straight-forward  relationship  between  IPS  and  the  spectral  contribution  of  any  single 
point  in  the  data  set  is  given  by  equation  (56).  Equation  (61)  defines  an  analogous  re¬ 
lationship  between  (57)  and  (59).  Unfortunately,  this  relationship  is  not  so  straight¬ 
forward  and  requires  additional  analysis  to  determine  the  benefits,  if  any,  to  be  gained 
from  processing  data  in  this  manner. 
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V.  RECOMMENDATIONS  AND  CONCLUSIONS 


The  comparative  behavior  of  the  spectrogram,  IPS,  WD  and  three  novel  t-f  distrib¬ 
utions  was  explored  in  Chapter  IV.  Of  the  three,  ImRD*  proves  to  be  particularly  sen¬ 
sitive  to  discontinuous  changes  in  frequency.  This  characteristic  behavior  may  be  useful 
in  certain  detection  applications.  When  confronted  with  more  complicated  spectra,  one 
containing  closely-spaced  stationary  components  or  continuously-varying  nonstationary 
components,  hnRDy  does  not  appear  particularly  useful.  Similarly  t##,!2  provides  an 
improvement  in  end-point  resolution  when  used  for  the  estimation  of  stationary  or 
pulsed  spectra.  Taking  the  square  root  of  |#/>j2  results  in  a  t-f  distribution,  denoted 
by  |  RD^f,!)  | ,  which  satisfies  many  of  the  desirable  properties  listed  in  Table  2.  In 
particular  |  /?£>(/»  |  : 

1.  Satisfies  both  zero  energy  requirements 

2.  Obeys  the  time  and  frequency  shift  properties 

3.  Is  positive  and  real  for  all  time  and  frequencies. 

How  well  i  RDf  1 2  can  track  a  rapidly  fluctuating  puised  signal,  similar  to  something 
found  in  frequency  hopped  communications,  is  an  area  worth  investigating. 

For  detection  of  continuously  changing  spectral  dy  namics.  !/?£>,!!  appears  to  be  a 
viable  processing  scheme,  however  the  performance  of  this  estimator  in  noise  needs  still 
to  be  examined.  Using  \RD, I2  as  an  estimator  is  not  without  problems.  Unlike 
I  RD\ ,  taking  the  square  root  of  does  not  produce  an  estimator  which  satisfies 

the  marginal  requit cments.  it  does  however  comply  with  the  zero  energy  and 
time  frequency  shift  properties  desirable  for  time-dependent  estimation  problems.  Used 
with  stationary  spectra  I RDt\l  appears  to  improve  the  resolution  of  closely  spaced 
components  beyond  that  currently  achieved  by  classical  methods.  For  stationary 
spectra,  1  RDji  is  a  biased  estimator,  i.s.,  the  true  spectral  peaks  occur  at  a  given  fixed 
offset  from  their  true  frequency  location,  PWD  also  provides  a  resolution  improvement; 
however  the  estimate  it  provides  never  settles  down  to  one  frequency  location.  For  short 
duration  data  \RDt\l  may  prove  superior  in  the  detection  of  multiple  stationary  com¬ 
ponents. 

In  addition  to  the  experimental  results  presented  in  Chapter  IV,  IPSy  demonstrates 
a  3-dB  noise  advantage  relative  to  PW'D  (Ref.  15J.  Coupled  with  superior  end-point 
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resolution  IPS  may  be  the  desirable  method  in  practical  analysis  problems.  All  the  re¬ 
sults  presented  were  derived  from  noise-free  data  sequences  using  digital  implementa¬ 
tions  (see  Appendix  A).  Because  ImRDy ,  |  RDy  \ 2  and  1  RDy  \  i  all  show  advantages  when 
applied  to  the  appropriate  signal  type,  analysis  of  their  noise  performance  is  an  open 
issue.  As  they  are  derived  from  the  RD,  as  is  IPS,  it  is  likely  that  they  enjoy  a  similar 
robustness. 

Looking  to  practical  applications  of  IPS,  a  brief  discussion  of  performance  in  a 
multi-sensor  environment  can  be  found  in  Appendix  B.  The  initial  results  look  promis¬ 
ing,  but  more  extensive  research  needs  to  be  conducted.  A  second  practical  application 
scheme  involves  the  use  of  a  cumulant  or  third-order  moment.  Typically  associated  with 
the  sonar  environment,  this  scheme  seeks  to  take  advantage  of  the  fact  that  the  odd 
moments  of  a  zero-mean,  gaussian  noise  process  are  identically  zero.  An  initial  investi¬ 
gation  using  IPS  to  compute  cumulants  can  be  found  in  Appendix  C. 

In  general,  WD  produces  very  narrow  spectral  ridges  but  suffers  form  poor  end¬ 
point  resolution  and  spectral  artifacts.  IPS  provides  an  improvement  over  these  short¬ 
comings  at  the  cost  of  spectral  broadening.  Recently  appearing  in  the  literature  is 
another  t-f  distribution  [Ref.  10,  22).  Defined  by  H.I.  Choi  and  W.J.  Williams,  this  t-f 
distribution  minimizes  the  effects  of  the  spectral  cross-terms.  Closer  examination  of  the 
resultant  spectra,  which  uses  the  kernel  function 

2  2 
7TL  T 

t)  =5  e  2a  , 

where  a  is  a  constant,  shows  that  the  spectral  ridges  are  broadened,  similar  in  behavior 
to  IPS.  In  classical  estimation,  i.e.,  the  periodogram,  one  maximizes  spectral  resolution 
using  nothing  but  the  raw  finite  data  set.  The  price  paid  is  a  large,  slow  roll-off  sidelobe 
structure  that  can  mask  a  true  component.  In  time-frequency  distributions,  i.e.,  in  the 
generalized  phase-space  equation  (19),  one  maximizes  spectral  detail  using  a  constant 
kernel  of  unit  amplitude.  The  price  paid  is  poor  end-point  resolution  and  spectral 
artifacts  across  time  and  frequency.  In  the  classical  analog,  using  any  other  window 
results  in  improved  sidelobe  behavior  at  a  corresponding  loss  in  detail  along  the  fre¬ 
quency  axis.  lPSy,  the  Rihaczek  derived  distributions  and  the  distribution  suggested  by 
Choi  and  Willirnas  all  improve  or  eliminate  the  disturbing  spectral  cross  terms  charac¬ 
teristic  of  WD.  The  price  paid  is  loss  in  spectral  detail. 
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APPENDIX  A.  COMPUTER  CODE 


To  conserve  space  only  the  read  file  and  four  of  the  si..  Fortran  codes  have  been 
included.  Data  was  generated  and  read  by  the  basic  programs  as  required.  The  data 
generation  file  is  not  included  here.  The  read  file  allows  easy  change  of  processing  pa¬ 
rameters.  Further,  IPSy  and  ImRDy  are  generated  by  code  that  difiers  in  a  minus  sign 
in  line  code.  lPSy  requires  a  plus  sign;  ImRDy  requires  a  minus  sign.  Similarly  |  1 2 

and  \RDy\l  are  simply  related.  The  location  of  the  sign  change  is  indicated  in  the  re¬ 
spective  algorithms.  Graphics  are  produced  using  DISSPLA. 

1.  Parameter  File 


1  1  031  039  08 

MODE  1 

HAMMING  WINDOW  ($' 

WTYPE 

B 

CONTR 

IPS: $' 

TTL 

FS=1/128 ,  40  POINTS  OF  DATA$ 1 

SIGNAL 

5  CELL  TIME  SMOOTHING? ' 

TSMTH 

BWLEN  EWLEN  WINC 
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OC^OOUOOOOOUOOOOOUOUUUOUOOOUUUUUUUUOOUCJOUOOUUUOOUOOOO 


2.  Spectrogram 


C**  THIS  FORTRAN  FILE  COMPUTES  THE  SPECTROGRAM  OF  *** 

C**  A  DATA  SEQUENCE  *** 

* 
* 

INPUT  DATA  SEQUENCE  IS  READ  USING  FILEDEF  4,  AS  THE  COMPLEX  * 
ARRAY  DATA(L).  * 

* 

L  IS  THE  LENGTH  OF  THE  DATA  SEQUENCE  AND  IS  ADJUSTED  FROM  THE  * 
PARAMETER  STATEMENT.  L  MUST  NOT  EXCEED  128.  * 

* 

ANALYSIS  PARAMETERS  ARE  READ  USING  FILEDEF  41.  THE  PARAMETERS  * 
ARE:  * 

ARGUMENT  TYPE  ALLOWED  VAULUES  * 

* 

MODE  -  II  1  PLOT  0  TO  PI  * 

2  PLOT  PI  TO  PI  * 

* 

PLTR  -  II  0  SHERPA  LASER  PRINTER 


1  IMB79  GRAPHICS  TERMINAL 


BWLEN-  13  3  DIGIT  INITIAL  WINDOW  LENGTH,  * 

MUST  BE  AN  ODD  INTEGER  * 

EWLEN-  13  3  DIGIT  FINAL  WINDOW  LENGTH  * 

WINC  -  12  2  DIGIT  WINDOW  INCREMENT,  MUST  * 

BE  AN  EVEN  INTEGER  * 

* 

WTYPE-  A19  19  CHARACTER  STRING  USED  IN  THE  * 

PLOT  HEADER  DISCRIBING  THE  * 

WINDOW  USED.  THE  CURRENT  * 

WINDOW  LENGTH  IS  AUTOMATICALLY  * 
INCLUDED  * 

* 

CONTR-  A1  1  CHARACTER  STRING  INDICATING  * 

TYPE  OF  PLOT  DESIRED  * 

* 

A  AMPLITUDE  PLOT  ONLY  * 

C  CONTOUR  PLOT  ONLY  * 

B  BOTH  AMPLITUDE  AND  CONTOUR  * 

* 

TTL  -  A43  43  CHARACTER  STRING  USED  IN  THE  * 

HEADING  WHICH  DESCRIBES  THE  * 

ALGORITHM  AND  THE  CLASS  OF  * 

SIGNAL  USED  * 

* 

SIGNAL-  A43  43  CHARACTER  STRING  DESCRIBING  * 

TEST  SIGNAL  * 

* 

vY 

OUT  REAL  * 

OUTPUT  ARRAY  OF  DIMENSION  512  BY  L  * 

COEF  COMPLEX 

ARRAY  OF  DATA  AFTER  IT  HAS  BEEN  WEIGHTED  WITH  A  * 

SLIDING  WINDOW  FUNCTION.  * 


UUOUUUUUUUUUUUUUUU  O  CJ  uooo 


* 

FT  COMPLEX  * 

ARRAY  OF  THE  1024  POINT  TRANSFORM  COEF  * 

* 

Z  INTEGER  * 

LENGTH  OF  CURRENT  WINDOW  * 

* 

M  INTEGER  * 

MID-POINT  OF  THE  CURRENT  WINDOW  * 

* 

AMAX  REAL  * 

MAXIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 

AMIN  REAL  * 

MINIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 


* 


PARAMETER(L=  32) 

INTEGER  I , J , N , M , MODE , Z , BWLEN , EWLEN , WINC , PLTR 
REAL  OUT(512,L) , AMAX, AMIN 
CHARACTER  WTYPE* 1 9 , TTL*43 , S I GNAL*43 , CONTR* 1 
COMPLEX  DATA( L) ,FT( 1024) ,COEF( 1024) 

CALL  EXCMSC'FILEDEF  4  DISK  TEST  IN  A  (PERM') 
CALL  EXCMSC'FILEDEF  41  DISK  PARAM  IN  A  (PERM') 


---  READ  IN  PARAMETER  LIST  . 

READ( 41,400) MODE , PLTR , BWLEN , EWLEN , WINC , WTYPE , CONTR , TTL , 

+  SIGNAL 

400  FORMAT  (1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43/ 
+  1X.A43) 


TEST  TO  ENSURE  WINDOW  LENGTH  IS  APPROPRIATE  . 

I=L-1 

IF  ((BWLEN  .GT.  I)  .OR.  (EWLEN  .  GT.  I))  THEN 
WRITE(*,69) 

GO  TO  99 
ENDIF 

I=M0D(BWLEN,2) 

K=MOD(WINC ,2) 

IF  (I  .EQ.  0)  THEN 

IF  (K  .EQ.  1)  THEN 
WRITE(*,68) 

GO  TO  99 

ELSE 

WRITE(*,67) 

GO  TO  99 
ENDIF 

ENDIF 

69  FORMAT  (IX, 'WINDOW  LENGTH  EXCEEDS  LENGTH  OF  THE  DATA') 

68  FORMAT  (IX,1 WINDOW  INCREMENT  MUST  BE  EVEN') 

67  FORMAT  (IX, 'INITIAL  WINDOW  LENGTH  MUST  BE  ODD') 

C . 
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o  o  or: 


C 

C . PLOTTING  DEVICE  CALL 

IF  (PT-TR  .EQ.  0)  THEN 
CALL  COMPRS 
ELSE 

CALL  IEM79 
END  IF 


PI=4*ATAN( 1.  ) 
READ(4,*)(D^TA(I) ,I=1,L) 


DO  111  Z=BWLEN , EWLEN , WINC 

M=(Z-l)/2 

CALL  ANGLE (0.  0) 

AMAX=0 
AMIN=AMAX 
DC  10  J=1,L 
DO  20  N=-M,M 

IF  (  ((I+N)  .GE.  1)  .AND.  ((I+N)  .  LE.  L)  )  THEN 
C0EF(M+N+1)=DATA( I+N) 

+  *(0. 54+0.  46*COS( 2*PI*N/( 2*M) ) ) 

ELSE 

COEF(M+N+l)=(C.  ,0.  ) 

END  IF 

20  CONTINUE 

DO  30  Nr-2*M+2,  J.024 
COEF(N)=(0.  ,0.  ) 

30  CONTINUE 

CALL  FFT(1024,C0EF,FT) 

IF  (  MODE.  EQ.  2)  THEN 
DO  40  N=l, 513,2 

OUf( INT((N+l)/2+255),I)=ABS(FT(N))**2/(2*M+l) 
IF  (OUT( INT( (N+J )/2+255) , I)  . GT.  AMAX)  THEM 
AMAX=0UT( INT( ( Nf 1 ' /2+255 ) , I ) 

ENLTF 

IF  (0  T( TNT((N+l)/2+255) ,1)  . LT.  AMIN)  THEN 
AMIN=0CT( INT( ( N+l) / 2+255 ) , I ) 

END  IF 

40  CONTINUE 

DO  50  N-515 , 1024,2 

OUT(INT((N-:)/2-256),I)=ABS(FT(N))**?/(2*M+l) 
IF  (OUT(INT((N-l)/2-256),l)  .  GT.  AMAX)  THEN 
AMAX=OUT( INT( ( N- 1 ) / 2 - 25  6 ) , I ) 

END  IF 

IF  (0UT( INT( (N  l)/2-256),I)  .  LT.  AMIN)  THEN 
AMIN=OUT(INT( ( N - 1 ) / 2 -2 5 6 ) , I ) 

ENDIF 

50  CONTINUE 

ELSE 

DO  51  N=l,512 

OUT(N,I)=ABS(FT(N',)**Z/(2*M+l) 

IF  (0UT(N,I)  .GT.  AMAX)  THEN 
AMAX=OUT(N,I) 

ENDIF 
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noon  oooo 


IF  (OUT(N, I)  .LT.  AMIN)  THEN 
AMIN=0UT(N, I) 

ENDIF 

51  CONTINUE 

ENDIF 

10  CONTINUE 

IF  ( CGNTR  .EQ.  'C')  THEN 
GO  TO  98 
ENDIF 

CALL  AREA2D(8.  ,9.  ) 

CALL  V0LM3D( 10.  ,10.  ,8.  ) 

CALL  HEADIN(TTL, 100,1. ,3) 

CALL  HEADIN( SIGNAL, 100,1. ,3) 

CALL  HEADINCNO  TIME  SMOOTHING,  1024  FFT$', 100,1.  ,3) 
CALL  MESSAGC WTYPE, 100, 2.  5,9.  3) 

CALL  INTN0(2*M+1, f ABUT' , 'ABUT1 ) 

CALL  MESSAGC  P0INTS)$' ,100, 'ABUT' ,' ABUT') 

CALL  X3NAME(' FREQUENCY  AXIS$’,100) 

CALL  Y3NAMF.C  'TIME  AXIS$',100) 

CALL  Z3NAME(’  $',100) 

CALL  VUANGL(-65.  ,70.  ,700. ) 

CALL  XNONUM 
C  CALL  ZNONUM 

CALL  MX1ALF( ' STANDARD' ,'#') 

CALL  MX2ALF( 'L/CGREEK' , '+' ) 

CALL  ANGLE( -25. 0) 

IF  (  MODE  .EQ.  1  )  THEN 
CALL  MESSAGC'  +0 #  ',5,0. ,2.3) 

ELSE 

CALL  MESSAGC'  +-P#  ',6,0. ,2.3) 

ENDIF 

CALL  ANGLE (-25.0) 

CALL  MESSAGC'  +P#  ’,5,4.9,0.15) 

CALL  GRAF3DC -256,256,256,0,L,L, 1.  0*AMIN, 

+  0. 5*( AMAX+AMIN) , 1.  0*AMAX) 

CALL  SURMAT(OUT,512,512,1,L,0.  ) 

CALL  ENDPL(O) 

C 

98  IF  ( CONTR  . NE.  'A')  THEN 

CALL  CONTORC  OUT , 5 1 2 , L , TTL , S I GNAL , WTYPE , Z , AMAX ) 


ENDIF 

111 

CONTINUE 

CALL  DONEPL 

99 

STOP 

END 

SSSSSSSSSSSFSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE  CONTORC  A , NX , NY , TITLE , S IGNL , VNDW , WLEN , AMAX ) 

THTS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS. 
NOTE:  THE  ARRAY  MUST  BE  REAL* 4. 

A  :  SINGLE  PRECISION  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS 
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NX:  NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NY:  NUMBER  OF  POINTS  IN  THE  Y-DIRECTION 
ZINC:  CONTOUR  INTERVAL 

DIMENSION  A(NX,NY) 

COMMON  WORK( 50000) 

SET  PARAMETERS  FOR  AXES: 

XORIG--256. 

XSTP=256. 

XMAX=256. 

Y0RIG=0. 

YSTP=NY 

YMAX=NY 

SET  CONTOUR  LEVEL 
ZINC=AMAX/10. 

CALL  SETCLR('CYAN') 

SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT: 

CALL  PAGEC8. 5,11. 0) 

CALL  BCOMON( 50000) 

CALL  AREA2D(6. 0,8. 0) 

LABEL  AXES: 

CALL  XNAME( 'FREQUENCY  -  AXIS  $',100) 

CALL  YNAME( 'TIME  -  AXIS  $’,100) 

CALL  GRAF( XORIG ,XSTP , XMAX , YORIG , YSTP , YMAX) 

CALL  FRAME 

TITLE: 

CALL  HEADINC CONTOUR  PLOT$' ,100,1.  ,3) 

CALL  HEADIN(TITLE, 100,1. ,3) 

CALL  HEADIN(SIGNL, 100,1. ,3) 

CALL  ANGLE(0.  0) 

CALL  MESSAG(WNDW, 100,1. 5,-. 7) 

CALL  INTNO(WLEH  ,' ABUT’ ,’ ABUT' ) 

CALL  MESSAG('  POINTS) $’ ,100, 'ABUT' ,' ABUT') 

MAKE  CONTOURS  AND  DRAW: 

CALL  SETCLR('RED') 

CALL  C0NMIN(3.  0) 

CALL  C0NANG( 60.  ) 

CALL  CONLIN(0,'MYCON’ ,' NOLABELS' ,2,10) 

CALL  CONMAK( A , NX , NY , ZINC ) 

CALL  CONTUR( 1 , ' LABELS ' , ’ DRAW ' ) 

CALL  ENDPL(O) 

RETURN 

END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 
SUBROUTINE  MYCON( RARAY , IARAY) 
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DIMENSION  RARAY ( 2 ) , I ARAY ( 1 ) 


THIS  ROUTINE  MAKES  NEGATIVE  CONTOURS  DASHED  AND  THE  ZERO  LINE  HEAVIER. 

CALL  RESET('DASH') 

IF  (RARAY(l)  .GE.  0. )  GO  TO  10 
CALL  DASH 
10  RARAY(2)  =  1. 

IARAY(l)  =  1 

IF  (RARAY(l)  .EQ.  0.)  IARAY(l)  =  2 

RETURN 

END 

JSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


****************************************************** 


*  * 

*  CALL  FFT(N,XTEMP,X)  * 

*  * 

*  X  -  OUTPUT  COMPLEX  ARRAY  CONTAINING  FFT  (1024)  * 

*  N  -  NUMBER  OF  POINTS  * 

*  XTMP  -  COMPLEX  ARRAY  CONTAINING  DATA  SAMPLES  * 

*  (starting  at  l,up  to  1024)  * 


****************************************************** 


SUBROUTINE  FFT(N,XTMP,X) 

COMPLEX  X( 1024) ,XTMP( 1024) ,WTFAC,TMP 
M=INT( L0G10(FL0AT(N) ) /L0G10( 2.  )+0.  5 ) 

EN  =  N 

PI  =  4.  0*ATAN(  1.0) 

DO  10  K=0 ,N-1 
NEWADR  =  0 
MADDR  =  K 
DO  20  1=0, M-l 

LRMNDR  =  MOD (MADDR, 2) 

NEWADR  =  NEWADR  +  LRMNDR*2**(M-1-I) 
MADDR  =  MADDR/ 2 
20  CONTINUE 

X(NEWADR+1)  =  XTMP(K+1) 

10  CONTINUE 
DO  50  L=1,M 

I SPACE  =  2**L 
S  =  N/ I SPACE 
IWIDTH  =  ISPACE/2 
DO  40  J=0 , ( IWIDTH-1) 

R  =  S*J 

ALPHA  =  2. *PI*R/EN 

WTFAC  =  CMPLX(  COS(ALPHA),  -SIN(ALPHA)) 
DO  30  ITOP=J ,N-2 , ISPACE 
IBOT  =  ITOP  +  IWIDTH 
TMP  =  X(IBOT+l)*WTFAC 
X( IBOT+1)  =  X( ITOP+1)  -  TMP 
X( ITOP+1)  =  X( ITOP+1)  +  TMP 
30  CONTINUE 

40  CONTINUE 
50  CONTINUE 


FFT00130 

FFT00140 

FFT00210 

FFT00270 

FFT00320 

FFT00330 

FFT00340 

FFT00350 

FFT00360 

FFT00370 

FFT00380 

FFT00390 

FFT00400 

FFT00410 

FFT00530 

FFT00610 

FFT00620 

FFT00630 

FFT00670 

FFT00720 

FFT00730 

FFT00740 

FFT00750 

FFT00S00 

FFT00810 

FFT00820 

FFT00830 

FFT00840 

FFT00850 

FFT00860 
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RETURN 

END 


FFT01000 

FFTOIOIO 
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3, 


IPS 


C**  THIS  FORTRAN  FILE  COMPUTES  THE  IPS  OF  A  DATA  SEQUENCE  *** 

* 

* 

INPUT  DATA  SEQUENCE  IS  READ  USING  FILEDEF  4,  AS  THE  COMPLEX  * 
ARRAY  X(N).  * 

v- 

N  IS  THE  LENGTH  OF  THE  DATA  SEQUENCE  AND  IS  ADJUSTED  FROM  THE  * 
PARAMETER  STATEMENT.  N  MUST  NOT  EXCEED  128.  * 

* 

ANALYSIS  PARAMETERS  ARE  READ  USING  FILEDEF  41.  THE  PARAMETERS  * 

ARE:  * 

ARGUMENT  TYPE  ALLOWED  VAULUES  * 

* 

MODE  “II  1  PLOT  0  TO  PI  * 

2  PLOT  PI  TO  PI  * 

* 

PLTR  -  II  0  SHERPA  LASER  PRINTER  * 

1  IMB79  GRAPHICS  TERMINAL  * 

'  * 

BWLEN-  13  3  DIGIT  INITIAL  WINDOW  LENGTH,  * 

MUST  BE  AN  ODD  INTEGER  * 

EWLEN-  13  3  DIGIT  FINAL  WINDOW  LENGTH  * 

WINC  -  12  2  DIGIT  WINDOW  INCREMENT,  MUST  * 

BE  AN  EVEN  INTEGER  * 

* 

WTYPE-  A19  19  CHARACTER  STRING  USED  IN  THE  * 

PLOT  HEADER  DISCRIBING  THE  * 

WINDOW  USED.  THE  CURRENT  * 

WINDOW  LENGTH  IS  AUTOMATICALLY  * 
INCLUDED  * 

* 

CONTR-  A1  1  CHARACTER  STRING  INDICATING  * 

TYPE  OF  PLOT  DESIRED  * 

it 

A  AMPLITUDE  PLOT  ONLY  * 

C  CONTOUR  PLOT  ONLY  * 

B  BOTH  AMPLITUDE  AND  CONTOUR  * 

* 

TTL  -  A43  43  CHARACTER  STRING  USED  IN  THE  * 

HEADING  WHICH  DESCRIBES  THE 
ALGORITHM  AND  THE  CLASS  OF  * 

SIGNAL  USED  * 

it 

SIGNAL-  A43  43  CHARACTER  STRING  DESCRIBING  * 

TEST  SIGNAL  * 

it 

TSMTH-  A25  25  CHARACTER  STRING  DESCRIBING  * 

TYPE  OF  TIME  SMOOTHING  USED  * 

it 
it 

OUT  REAL 

OUTPUT  ARRAY  OF  DIMENSION  M2  BY  N  * 

C  * 

C  SAMp  COMPLEX  * 
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SHIFTED  VERSION  OF  X 
SAM  COMPLEX 

SHIFTED  AND  CONJUGATED  VERSION  OF  X 
C  COMPLEX 

ARRAY  OF  SUM  OF  PRODUCTS  OF  DIMENSION  1  BY  1025 
POSITIONS  1  -  512  ARE  CONJUGATE  SYMMETRIC  WITH 
POSITIONS  514  -  1025.  POSITION  513  =  (0. ,0. ) 

FT  COMPLEX 

ARRAY  OF  THE  1024  POINT  TRANSFORM  C 

Z  INTEGER 

LENGTH  OF  CURRENT  WINDOW 

M  INTEGER 

MID-POINT  OF  THE  CURRENT  WINDOW 
AMAX  REAL 

MAXIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS 
AMIN  REAL 

MINIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS 

NX  HORIZONTAL  DIMENSION  OF  OUT  WHICH  IS  ALWAYS  512 

DATA  CAN  BE  OUTPUT  USING  FILEDEF  61.  POINTS  OF  INTEREST  MUST 
BE  DEFINED  IN  THE  APPROPRAITE  SECTION  OF  CODE.  BECAUSE 
OF  SPACE  CONSTRAINTS,  THE  DATA  OUTPUT  FILE  IS  WRITTEN  TO 
THE  B  DISK. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

it 


^r*************************************- 


PARAMETERS  64) 

COMPLEX  X(N) ,C( 1025) , SAMP, SAM, FT( 1024) 

REAL  0UT(512,N) , AMAX, AMIN 

INTEGER  NX , K , I , J , MODE ,Z,M, BWLEN , EWLEN , WINC , PLTR 
CHARACTER  WTYPE*19 ,TTL*43 ,SIGNAL*43 ,TSMTH*25 ,CONTR*l 
CALL  EXCMS( 'FILEDEF  4  DISK  TEST  IN  (PERM') 

CALL  EXCMS( 'FILEDEF  41  DISK  PARAM  IN  (PERM') 

CALL  EXCMS( 'FILEDEF  61  DISK  DATA  OUT  B  (PERM') 

.  READ  IN  PARAMETER  LIST  . 

RE AD( 4 1 , 400 ) MODE , PLTR , BWLEN , EWLEN , WINC , WTYPE , CONTR , TTL , 

+  SJGNAL,TSMTH 

FORMAT  (1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X, A43 
+  /IX,A43/1X,A25) 


....  JEST  TO  ENSURE  WINDOW  LENGTH  IS  APPROPRIATE 
I=N-1 

IF  ((BWLEN  .GT.  I)  .OR.  (EWLEN  . GT. I))  THEN 

«tnr  Tr  f  k  n  \ 

1L(  -  ,  uy  ) 

GO  TO  99 
END  IF 
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I=M0D(BWLEN,2) 

K=M0D(WINC,2) 

IF  (I  .EQ.  0)  THEN 

IF  (K  .EQ.  1)  THEN 
WRITE(*,68) 

GO  TO  99 

ELSE 

WRITE(*,67) 

GO  TO  99 
END  IF 

END  IF 

69  FORMAT  (IX,' WINDOW  LENGTH  EXCEEDS  LENGTH  OF  THE  DATA') 

68  FORMAT  (IX,' WINDOW  INCREMENT  MUST  BE  EVEN') 

67  FORMAT  (IX, 'INITIAL  WINDOW  LENGTH  MUST  BE  ODD') 


---PLOTTING  DEVICE  CALL 
IF  (PLTR.EQ. 0)THEN 
CALL  COMPRS 

ELSE 

CALL  IBM79 
ENDIF 


PI=4*ATAN( 1.  ) 

NX=512 

READ  (4,*)(X(J),J=1,N) 

DO  111  Z=  BWLEN , EWLEN , WINC 
WRITE( 6 1 , 600 )TTL , S IGNAL , TSMTH , WTYPE , Z 
600  FORMAT  ( 1X,A43/1X,A43/1X,A25/1X,A19,I3, '  POINTS)') 

CALL  ANGLE(0.  ,0.  ) 

M=( Z - 1 ) / 2 
AMAX=0. 

AMIN=AMAX 
DO  10  1=1. N 
DO  20  K=0,512 
SAMP=(0. ,0. ) 

SAM=SAMP 

IF  (  ( I+K)  .  LE.  N  )  THEN 
SAMP=X( I+K) 

ENDIF 

IF  (  ( I-K)  .GT.  0  )  THEN 
SAM=CONJG(X( I-K) ) 

ENDIF 

C 

C+  -  +  -  +  “  +  -  +  -  +  -  +  -  +  “  +  -  +  “  +  -  +  -  +  "  +  -  +  -  +  “  +  -  + 
C+-+-+-+  sum  of  product  is  IPS  difference  of  product  is  IMRD  -+-+-+ 
IF  (K  .LE.  M)  THEN 

C(K+1)=(X( I)*SAM  +  CONJG(X( I ) )*SAMP) 

+  *(0.  54+0.  46*C0.'i<'  2*PI*K/(  2*M) ) ) 

ELSE 

C(K+1)=0 

hNUil 

C+  -+-+-+-+-+-+-+-+-+-+-+-+-  +-+-+-+-+ 
C 


81 


nnn  o  o 


C(1024-K+1)=CONJG(C(K+1)) 

20  CONTINUE 

C(513)=(0.  ,0.  ) 

CALL  FFT( 1024,  :,FT) 

DO  40  K=1 ,M0DE*j 12 ,M0DE 
IF  (REAL(FT(K))  .  GT.  AMAX)  THEN 
AMAX=REAL(FT(K) ) 

END  IF 

IF  (REAL(FT(K))  . LT.  AMIN)  THEN 
AMIN=REAL(FT(K)) 

ENDIF 

IF  (  (K  .LT;  514)  .AND.  (MODE  .  EQ.  2))  THEN 
OUT( INT( (K+l)/ 2+255) , I)=REAL(FT(K) ) 

ELSE 

IF  (  MODE  .EQ.  2  )  THEN 

OUT( INT( ( K+l ) / 2 -25 7 ) , I )=REAL( FT( K) ) 
ELSE 

OUT( K , I )=REAL( FT( K) ) 

ENDIF 

ENDIF 


40  CONTINUE 

10  CONTINUE 

C . FOR  TIME  SMOOTHING  PURPOSES  . 

DO  48  K=l,512 
DO  46  1=1, N-2 

OUr(K,I)=(OUT(K,I)+OUT(K,I+l)+OUT(K  1+2) )/3 

46  CONTINUE 

DO  47  I=N,3,-1 

OUT(K,I)=(OUT(K,I) +OUT ( K , I - 1 ) +OUT (K, 1-2) )/3 

47  CONTINUE 
C 

C . DATA  OUTPUT . 

IF((K.  GE.  208).  AND. (K. LT. 370))THEN 
C  IF( ( (K.  GE. 120). AND. (K. LT. 140)). OR. ((K. GE. 370). AND. 

C  +  (K.  LT.  390)))THEN 

WRITE(61,601) 

DO  81  I=1,N-1, 14 

WRITE(61,602)K,I,0UT(K,I) 

81  CONTINUE 

ENDIF 

601  FORMAT  ('FREQ  BIN=' ,8X, ’TIME  BIN=’ ,7X, 'AMPLITUDE^ ) 

602  FORMAT  ( 10X, 14, 13X, 14, 14X,E14.  7) 


48 


CONTINUE 


.  PLOTTING 

IF(CONTR. EQ.  *  C  * )THEN 
GO  TO  50 
ENDIF 

CALL  BSHIFT  (  -0.2  ,-.25) 
CALL  AREA2D( 8,9) 

CALL  VOLM3DC 10,10,8) 

CALL  HEADIN(TTL, 100,1.  ,3) 
CALL  HEADIN( SIGNAL,  100,1..  ,3) 
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CALL  HEADIN(TSMTH, 100,1.  ,3) 

CALL  MESSAG( WTYPE , 100,2.  5 }9.  3) 

CALL  INTN0(Z, 'ABUT'  ’ABUT1) 

CALLMESSAGC  POINT..  )$', 100 .' ABUT V  ABUT' ) 

CALL  X3NAMEC FREQUENCY  AXIS$f,100) 

CALL  Y3NAME( ' TIME  AXIS$',100) 

CALL  Z3NAME( '  $’,100) 

CALL  VUANGL(-65,70,700) 

CALL  XNONUM 
C  CALL  ZNONUM 

CALL  MX1ALF(  ’  STANDARD'  ,'#') 

CALL  MX2ALF( '  L/CGREEK' , '  •*•' ) 

CALL  ANGLE(-25.0) 

IF  (  MODE  .EQ.  2  )  THEN 

CALLMESSAGC  +-P#  ',6,0.  ,2.3) 

ELSE 

CALLMESSAGC  +0#  ',5,0.  ,2.  3) 

ENDIF 

CALL  ANGLEC-25.0) 

CALL  MESSAGC  +P #  ’,5,4.9,0.15) 

CALL  GRAF3D( -256,256,256, 1 , N , N , 1 .  0*AMIN , 

+  0. 5*(AMAX-AMIN),1.0*AMAX) 

CALL  SURMATC OUT ,522, 512,1, N,0. ) 

CALL  ENDPL(O) 

C 

50  IF  (CONTR.NE. 'A')THEN 

DO  49  1=1, N 

DO  51  K=l,512 

IF  (OUT(K,I)  .LT.  0)  THEN 
OUT(K,I)=0 
ENDIF 

51  CONTINUE 

49  CONTINUE 

CALL  CONTOR ( OUT , NX , N , TTL , S I GNAL , WTYPE , TSMTH , Z , AMAX ) 

ENDIF 

C . 

c 

WRITE( 6 1,603) AMAX, AMIN 
603  FORMAT  (IX, 'MAXIMUM  AMPLITUDE=' ,E14. 7 
+  /'MINIMUM  AMPLITUDE^ ,E14.  7) 

111  CONTINUE 

CALL  DONEPL 
99  STOP 
END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE  CONTOR (A , NX , NY , TITLE .SIGNL, WNDW , TAVG , WLEN , AMAX) 

THIS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS. 
NO  IT:  THE  ARRAY  MUST  BE  REAL*4. 

A  :  SINGLE  PRECISION  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS 
NX:  NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NY:  NUMBER  OF  POINTS  IN  THE  Y-DIRECTION 
ZINC:  CONTOUR  INTERVAL 
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DIMENSION  A(NX,NY) 

COMMON  W0RK( 50000) 

C 

C  SET  PARAMETERS  FOR  AXES: 

X0RIG=-256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

C 

C  SET  CONTOUR  LEVEL 
ZINC=AMAX/10. 

C 

CALL  SETCLRC ’CYAN' ) 

C 

C  SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT: 

CALL  PAGE( 8.  5,11.0) 

CALL  BCOMON(SOOOO) 

CALL  AREA2D(6.  0,8.  0) 

LABEL  AXES: 

CALL  XNAME( 'FREQUENCY  -  AXIS  $',100) 

CALL  YNAMEC'TIME  -  AXIS  $’,100) 

CALL  GRAF( XORIG , XSTP , XMAX , YORIG , YSTP , YMAX) 

CALL  FRAME 

TITLE: 

CALL  HEADINC CONTOUR  PL0T$' ,100,1.  ,4) 

CALL  HEADINC TITLE, 100,1. ,4) 

CALL  HEADINC SIGNL, 100,1.  ,4) 

CALL  HEADINCTAVG, 100,1. ,4) 

CALL  ANGLECO.  0) 

CALL  MESSAGCWNDW, 100,1. 5 , -.  7) 

CALL  INTNOCWLEN  , 'ABUT' , 'ABUT' ) 

CALL  MESSAGC'  POINTS) $' ,100, 'ABUT' ABUT') 

MAKE  CONTOURS  AND  DRAW: 

CALL  SETCLRC 'RED') 

CALL  C0NMINC3. 0) 

CALL  C0NANGC60. ) 

CALL  CONLINCO.'MYCON' ,’ NOLABELS' ,2,10) 

CALL  CONMAKC A , NX , NY , ZINC ) 

CALL  CONTURC 1 , ' LABELS ' , ' DRAW ’ ) 

CALL  ENDPLCO) 

RETURN 
END 

ISSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE  MYCON ( RARAY , I ARAY ) 

DIMENSION  RARAYC2) ,IARAY( 1) 

THIS  ROUTINE  MAK.  NEGATIVE  CONTOURS  DASHED  AND  THE  ZERO  LINE  HEAVIER. 
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CALL  RESET('DASH’) 

IF  (RARAY(l)  .GE.  0.)  GO  TO  10 
CALL  DASH 
10  RARAY(2)  =  1. 

IARAY(l)  =  1 

IF  (RARAY(l)  .EQ.  0.)  IARAY(l)  =  2 
RETURN 
END 
C 

csssssssssssssssssssssssssssssssssss 

*5 ***************************************************** 


A 

A 

*  CALL  FFT(N,XTEMP,X) 

A 

A 

*  X  -  OUTPUT  COMPLEX  ARRAY  CONTAINING  FFT  (1024) 

A 

*  N  -  NUMBER  OF  POINTS 

A 

*  XTMP  -  COMPLEX  ARRAY  CONTAINING  DATA  SAMPLES 

A 

*  (starting  at  l,up  to  1024) 

A 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SUBROUTINE  FFT(N,XTMP,X) 

FFT00130 

COMPLEX  X( 1024) ,XTMP( 1024) ,WTFAC,TMP 

FFT00140 

M=INT(L0G10(FL0AT(N) )/LOG10( 2. )+0. 5) 

EN  =  N 

FFT00210 

PI  =  4. 0*ATAN( 1. 0) 

FFT0027 

DO  10  K=0 ,N- 1 

FFT00320 

NEWADR  =  0 

FFT00330 

MADDR  =  K 

FFT00340 

DO  20  1=0, M-l 

FFT00350 

LRMNDR  =  MOD (MADDR, 2) 

FFT00360 

NEWADR  =  NEWADR  +  LRMNDR*2**(M-1-I) 

FFT00370 

MADDR  =  MADDR/ 2 

FFT00380 

20  CONTINUE 

FFT00390 

X(NEWADR+1)  =  XTMP(K+1) 

FFT00400 

10  CONTINUE 

FFT00410 

DO  50  L=1,M 

FFT00530 

I SPACE  =  2**L 

FFT00610 

S  =  N/ 1  SPACE 

FFT00620 

IWIDTH  =  ISPACE/2 

FFT00630 

DO  40  J=0,( IWIDTH-1) 

FFT00670 

R  =  S*J 

FFT00720 

ALPHA  =  2. *PI*R/EN 

FFT00730 

WTFAC  =  CMPLX(  COS( ALPHA),  -SIN(ALPHA)) 

FFT00740 

DO  30  ITOP=J , N- 2 , ISPACE 

FFT00750 

IBOT  =  ITOP  +  IWIDTH 

FFT00800 

TMP  =  X(IBOT+l)*WTFAC 

FFT00810 

X( IBOT+1)  =  X( ITOP+1)  -  TMP 

FFT00820 

X( ITOP+1)  =  X( ITOP+1)  +  TMP 

FFT00830 

30  CONTINUE 

FFT00840 

40  CONTINUE 

FFT00850 

50  CONTINUE 

FFT00860 

RETURN 

FFT01000 

END 

FFT01010 
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CVc* 

THIS  FORTRAN  FILE  COMPUTES 

i  THE  PWD  OF  A  DATA  SEQUENCE  *** 

C 

* 

C 

* 

C 

INPUT  DATA  SEQUENCE  IS  READ  USING  FILEDEF  4,  AS  THE  COMPLEX 

* 

c 

ARRAY  X(N). 

* 

c 

* 

c 

N  IS  THE  LENGTH  OF  THE  DATA  SEQUENCE  AND  IS  ADJUSTED  FROM  THE 

* 

c 

PARAMETER  STATEMENT.  N 

MUST  NOT  EXCEED  128. 

* 

c 

* 

c 

ANALYSIS  PARAMETERS  ARE  READ 

USING  FILEDEF  41.  THE  PARAMETERS 

* 

c 

ARE: 

* 

c 

ARGUMENT  TYPE 

ALLOWED  VAULUES 

* 

c 

* 

c 

MODE  -  11 

1  PLOT  0  TO  PI 

* 

c 

2  PLOT  PI  TO  PI 

it 

c 

•it 

c 

PLTR  -  11 

0  SHERPA  LASER  PRINTER 

it 

c 

1  IMB79  GRAPHICS  TERMINAL 

* 

c 

* 

c 

BWLEN-  13 

3  DIGIT  INITIAL  WINDOW  LENGTH, 

* 

c 

MUST  BE  AN  ODD  INTEGER 

it 

c 

EWLEN-  13 

3  DIGIT  FINAL  WINDOW  LENGTH 

* 

c 

WINC  -  12 

2  DIGIT  WINDOW  INCREMENT,  MUST 

K 

c 

BE  AN  EVEN  INTEGER 

it 

c 

it 

c 

WTYPE-  A 19 

19  CHARACTER  STRING  USED  IN  THE 

it 

c 

PLOT  HEADER  DISCRIBING  THE 

it 

c 

WINDOW  USED.  THE  CURRENT 

it 

c 

WINDOW  LENGTH  IS  AUTOMATICALLY 

it 

c 

INCLUDED 

it 

c 

it 

c 

CONTR-  A1 

1  CHARACTER  STRING  INDICATING 

it 

c 

TYPE  OF  PLOT  DESIRED 

it 

c 

it 

c 

A  AMPLITUDE  PLOT  ONLY 

it 

c 

C  CONTOUR  PLOT  ONLY 

it 

c 

B  BOTH  AMPLITUDE  AND  CONTOUR 

it 

c 

it 

c 

TTL  -  A43 

43  CHARACTER  STRING  USED  IN  THE 

it 

c 

HEADING  WHICH  DESCRIBES  THE 

it 

c 

ALGORITHM  AND  THE  CLASS  OF 

it 

c 

SIGNAL  USED 

it 

c 

it 

c 

SIGNAL-  A43 

43  CHARACTER  STRING  DESCRIBING 

it 

c 

p 

TEST  SIGNAL 

it 

«'f 

c 

TSMTH-  A25 

25  CHARACTER  STRING  DESCRIBING 

it 

c 

TYPE  OF  TIME  SMOOTHING  USED 

it 

c 

it 

c 

it 

c 

OUT  REAL 

it 

c 

p 

OUTPUT  ARRAY  OF  DIMENSION  512  BY  N 

it 

it 

c 

SAMP  COMPLEX 

it 
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SHIFTED  VERSION  OF  %  * 

* 

SAM  COMPLEX  * 

SHIFTED  AND  CONJUGATED  VERSION  OF  X  * 

* 

C  COMPLEX  * 

ARRAY  OF  SUM  OF  PRODUCTS  OF  DIMENSION  1  BY  1025  * 

POSITIONS  1-512  ARE  CONJUGATE  SYMMETRIC  WITH  * 

POSITIONS  514  -  1025.  POSITION  513  =  (0. ,0. )  * 

* 

FT  COMPLEX  * 

ARRAY  OF  THE  1024  POINT  TRANSFORM  C  * 

* 

Z  INTEGER  * 

LENGTH  OF  CURRENT  WINDOW  * 

* 

M  INTEGER  * 

MID-POINT  OF  THE  CURRENT  WINDOW  * 

* 

AMAX  REAL  * 

MAXIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 

AMIN  REAL  * 

MINIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 

NX  HORIZONTAL  DIMENSION  OF  OUT  WHICH  IS  ALWAYS  512  * 

* 


DATA  CAN  BE  OUTPUT  USING  FILEDEF  61.  POINTS  OF  INTEREST  MUST  * 
BE  DEFINED  IN  THE  APPROPRAITE  SECTION  OF  CODE.  BECAUSE  * 
OF  SPACE  CONSTRAINTS,  THE  DATA  OUTPUT  FILE  IS  WRITTEN  TO  * 
THE  B  DISK.  * 


>V******Vf** 


********** rt*****rtir**rt***V«V**** ******************** ********* 


PARAMETER(N=  64) 

COMPLEX  X( 512) ,C( 1025) ,SAMP,SAM,FT( 1024) 

REAL  OUT(512,N), AMAX, AMIN 

INTEGER  NX , K , I , J .MODE ,Z,M, BWLEN , EWLEN , WINC , PLTR 
CHARACTER  WTYPE*19,TTL*43,SIGNAL*43,TSMTH*25 .G0NTR*1 
CALL  EXCMS( 'FILEDEF  4  DISK  TEST  IN  (PERM') 

CALL  EXCMS( 'FILEDEF  41  DISK  PARAM  IN  (PERM') 

CALL  EXCMS( 'FILEDEF  61  DISK  DATA  OUT  B  (PERM') 

■ . .  READ  IN  THE  PARAMETER  LIST . 

READ( 4 1 , 400 ) MODE , PLTR , BWLEN , EWLEN , WINC , WTYPE , CONTR , TTL , 

+  SIGNAL, TSMTH 

400  FORMAT( 1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43 
+  /1X,A43/1X,A25) 


.  TEST  TO  ENSURE  WINDOW  LENGTH  IS  APPROPRIATE 

I=N-1 

IF  ((BWLEN  .GT.  I)  .OR.  (EWLEN  .  GT.  I))  THEN 
WRITE(*,69) 

GO  TO% 99 
ENDIF 
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I*M0D(BWLEN.2T 

K-M0D(WINC,2) 

IF  (I  .EQ.  0)  THEN 

IF  (K  .EQ.  1)  THEN 
WRITE(*t68) 

GO  TO  99 

ELSE 

WRITER, 67) 

GO  TO  99 
ENDIF 

ENDIF 

69  FORMAT  f IX, ’WINDOW  LENGTH  EXCEEDS  LENGTH  OF  THE  DATA’) 

68  FORMAT  (IX,* WINDOW  INCREMENT  MUST  BE  EVEN’) 

67  FORMAT  (IX,* INITIAL  WINDOW  LENGTH  MUST  BE  ODD') 


■ . PLOTTING  DEVICE  CALL 

IFO-LTR.  EQ.  0)THEN 
CALL  COMPRS 
ELSE 

CALL  IBM79 
ENDIF 


PI=4*ATAN( 1. ) 

NX=512 

READ(4,*)(X( J) ,J=1,2*N,2) 

• . DATA  INTERPOLATION 

DO  5  J=2 , 2*N , 2 
X(J)=(0.  ,0.) 

5  CONTINUE 

CALL  FFT(2*N,X,FT) 

DO  10  J=N/2+2 , 2*N -N/ 2+1 
FT(J)=(0.  ,0.  ) 

10  CONTINUE 

DO  20  J=1,2*N 

FT( J)=C0N JG( FT( J) ) 

20  CONTINUE 

CALL  FFT(2*N,FT,X) 

DO  30  J=1,2*N 

X(J)=CONJG(X(J))/N 
30  CONTINUE 


DO  111  Z=BWLEN , EWLEN , WINC 
WRITE( 61 ,600)TTL, SIGNAL, TSMTH,WTYP£,Z 
600  FORMAT( 1X,A43/1X,A37/1X,A25/1X,A19 ,13, 1  POINTS)’) 
M=(Z-l)/2 
AMAX=0 

AWTtl-AMAV 

DO  40  I=1,2*N,2 
DO  50  K=0,512 
SAMP=(0.  ,0.  ) 
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SAM=(0.  ,0.  ) 

IF  (  (I+K)  .LE.  2*N  )  THEN 
SAMP=X( I+K) 

END  IF 

IF  (  (I-K)  .GT.  0  )  THEN 
SAM=CONJG(X(I-K)) 

END  IF 


IF  (K  .LE.  2*M)  THEN 

C(K+1)=SAMP*SAM*(0.  54+0.  46*C0S(2*PI*K/(4*M))) 
ELSE 

C(K+1)=0 
END  IF 


C( 1024-K+1)=C0NJG(C(K+1) ) 

CONTINUE 
C(513)=(0.  ,0.  ) 

CALL  FFT( 1024, C, FT) 

DO  60  K= 1 , M0DE*5 1 2 , MODE 

IF  (REAL(FT(K))  .  GT.  AMAX)  THEN 
AMAX=REAL(FT(K) ) 

ENDIF 

IF  (REAL(FT(K))  .  LT.  AMIN)  THEN 
AMIN=REAL(FT(K)) 

ENDIF 

IF  (  (K.  LT.  514)  .AND.  (MODE  .EQ.  2))  THEN 
OUT( INT( ( K+ 1 ) / 2+255 ) , ( 1+1 ) / 2 ) =REAL( FT( K) ) 
ELSE 

IF(  MODE  .EQ.  2  )  THEN 
0UT(INT((K+l)/2-257),(I+l)/2)=REAL(FT(K)) 
ELSE 

0UT(K, ( I+1)/2)=REAL(FT(K) ) 

ENDIF 
ENDIF 
CONTINUE 
CONTINUE 

--FOR  TIME  SMOOTHING  PURPOSES  . 

DO  48  K=l,512 
DO  46  1=1, N-2 

OUT(K,I)=(OUT(K,I)+OUT(K,I+l)+OUT(K,I+2))/3 
CONTINUE 
DO  47  I=N ,3,-1 

0UT(K,I)=(0UT(K,I)+0UT(K,I-l)+0UT(K,I-2))/3 
CONTINUE 

.  DATA  OUTPUT  . 

IF((K.  GT.  285).  AND.  (K.  LT.  310))THEN 
WRITE(61,601) 

DO  81  1=1, N 

WRITE(6l,602)K,I,OUT(K,I) 

81  CONTINUE 
ENDIF 

601  FORMAT( 'FREQ  BIN=' ,8X,' TIME  BIN=' ,7X, 'AMPLITIDE=' ) 
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602  F0RMAT( 10X , 14 , 13X , 14 , 14X , E 14.  7) 
C . 

c 

48  CONTINUE 


.  PLOTTING  . 

IF  (CONTR.EQ.  'C')THEN 
GO  TO  52 
ENDIF 

CALL  BSHIFT  (  -0. 2  ,  -0. 25) 

CALL  AREA2D(8,9) 

CALL  V0LM3D( 10,10,8) 

CALL  HEADIN(TTL, 100,1. ,3) 

CALL  HEADINC SIGNAL, 100,1.  ,3) 

CALL  HEADINC TSMTH, 100,1. ,3) 

CALL  ANGLECO.  ,0.  ) 

CALL  MESSAG(WTYPE, 100.2. 5,9. 3) 

CALL  INTNO(Z  , 'ABUT' , ' ABUT' ) 

CALLMESSAGC’  POINTS  )$' ,100.' ABUT' ,' ABUT’) 

CALL  X3NAMEC 'FREQUENCY  AXIS$ ' , 100) 

CALL  Y3NAMEC 'TIME  AXIS$’,100) 

CALL  Z3NAMEC'  $’,100) 

CALL  VUANGLC -65,70; 700) 

CALL  XNONUM 
C  CALL  ZNONUM 

CALL  MX1ALFC ' STANDARD' , ' #' ) 

CALL  MX2ALFC ' L/CGREEK' , ' +' ) 

CALL  ANGLE( -25. 0) 

IF  (  MODE  .EQ.  2  )  THEN 
CALL  MESSAGC  +-P#  ',6,0.  ,2.  3) 

ELSE 

CALLMESSAGC'  +0#  ',5,0. ,2. 3) 

ENDIF 

CALL  ANGLE (-25.0) 

CALLMESSAGC’  +P #  *,5,4.9,0.15) 

CALL  GRAF3D(-256,256,256,1,N,N,1.0*AMIN,.5*CAMAX-AMIN), 

+1. 0*AMAX) 

CALL  SURMAT(0UT,512,512, 1,N,0. ) 

CALL  ENDPL(O) 

C 

52  IFCCONTR.NE. 'A’)THEN 
DO  49  1=1, N 

DO  51  K=l,512 

IF((OUT(K,I).  LT.  5. 0). AND. (OUT(K,I)  .  GT.  -5.0))THEN 
OUT(K, I)=0 
ENDIF 

51  CONTINUE 

49  CONTINUE 

CALL  CONTORC  OUT , NX , N , TTL , S I GNAL , WTYPE , TSMTH , Z ) 

ENDIF 


C 

WRITE(61,603)AMAX,AMIN 
603  FORMAT ( IX, 'MAXIMUM  AMPLITUDE=' ,E14.  7 
+  /IX, 'MINIMUM  AMPLITUDE=' ,E14.  7) 
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111  CONTINUE 
CALL  DONEPL 
99  STOP 
END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE  C0NT0R( A , NX , NY , TITLE , S IGNL , WNDW , TAVG , WLEN) 

THIS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS. 
NOTE:  THE  ARRAY  MUST  BE  REAL*4. 

A  :  SINGLE  PRECISION  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS 
NX:  NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NY:  NUMBER  OF  POINTS  IN  THE  Y-DIRECTION 
ZINC:  CONTOUR  INTERVAL 

DIMENSION  A(NX,NY) 

COMMON  WORK( 50000) 

C 

C  SET  PARAMETERS  FOR  AXES: 

XORIG=-256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

C 

C  SET  CONTOUR  LEVEL 
ZINC=AMAX/10 
C 

CALL  SETCLR('CYAN') 

C 

C  SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT: 

CALL  PAGEf  8.  5,11.0) 

CALL  BCOMON( 50000) 

CALL  AREA2D(6.  0,8.  0) 

C 

C  LABEL  AXES: 

CALL  XNAMEC FREQUENCY  -  AXIS  $’,100) 

CALL  YNAME('TIME  -  AXIS  $’,100) 

CALL  XNONUM 
C 

*  CALL  GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

CALL  FRAME 
C 

C  TITLE: 

CALL  HEADINC CONTOUR  PLOT$' ,100,1.  ,4) 

CALL  HEADINC TITLE, 100,1.  ,4) 

CALL  HEADINC SIGNL, 100,1.  ,4) 

CALL  HEADINC TAVG, 100,1.  ,4) 

CALL  ANGLE (0.0) 

CALL  MESSAGC WNDW, 100, 1.5.-.  7) 

Urtiili  VVIjCjIN  ,  ADVA  ,  tiDVl  ) 

CALL  MESSAGC ' POINTS) $ ' , 100 , ' ABUT' , ’ ABUT' ) 

C 

C  MAKE  CONTOURS  AND  DRAW: 
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CALL  SETCLR('RED') 

CALL  C0NMIN(3.0) 

CALL  C0NANG(60.  ) 

CALL  CONLIN(0,'MYCON’ ,' NOLABELS' ,2,10) 

CALL  CONMAK(A,NX,NY,' SCALE') 

CALL  C0NTUR( 1 , ' LABELS ' , ' DRAW ' ) 

CALL  ENDPL(O) 

RETURN 

END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE  MYCON( RARAY , IARAY) 

DIMENSION  RARAY(2),IARAY(1) 


THIS  ROUTINE  MAKES  NEGATIVE  CONTOURS  DASHED  AND  THE  ZERO  LINE  HEAVIER. 

CALL  RESET('DASH’) 

IF  (RARAY(l)  .GE.  0.)  GO  TO  10 
CALL  DASH 
10  RARAY(2)  ■  1. 

IARAY(l)  =  1 

IF  (RARAY(l)  .EQ.  C. )  IARAY(l)  =  2 

RETURN 

END 

ssssssssssssssssssssssssssssssssss 


A 

* 

A 


CALL  FFT(N,XTEMP,X) 


**************************** 

* 

* 

* 


*  X  -  OUTPUT  COMPLEX  ARRAY  CONTAINING  FFT  (1024)  * 

*  N  -  NUMBER  OF  POINTS  * 

*  XTMP  -  COMPLEX  ARRAY  CONTAINING  DATA  SAMPLES  * 

*  (starting  at  l,up  to  1024)  * 


SUBROUTINE  FFT(N,XTMP,X)  FFT00130 

COMPLEX  X( 1024) ,XTMP( 1024) ,WTFAC ,TMP  FFT00140 

M=INT(LOG10(FLOAT(N))/LOG10(2,  )+0.  5) 

EN  =  N  FFT00210 

PI  =  4.  0*ATAN( 1. 0)  FFT00270 

DO  10  K=0,N-1  FFT00320 

NEWADR  =  0  FFT00330 

MADDR  =  K  FFTQG340 

DO  20  1=0, M-l  FFT00350 

LRMNDR  =  MOD(MADDR,2)  FFT00360 

NEWADR  =  NEWADR  +  LRMNDR*2**(M-1-I)  FFT00370 

MADDR  =  MADDR/ 2  FFT00380 

20  CONTINUE  FFT00390 

X(NEWADR+1)  =  XTMP(K+1)  FFT00400 

10  CONTINUE  •  ‘  FFT00410 

DO  50  L=1,M  FFT00530 

ISPACE  =  2**L  FFT00610 
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S  =  N/ I SPACE  FFT00620 

IWIDTH  ■  ISPACE/2  FFT00630 

DO  40  J=0,( IWIDTH- 1)  FFT00670 

R  =  S*J  FFT00720 

ALPHA  =  2.  *PI*R/EN  FFT00730 

WTFAC  =  CMPLXC  COS( ALPHA),  -SIN( ALPHA))  FFT00740 

DO  30  ITOP=J,N-2,ISPACE  FFT00750 

I BOX  =  ITOP  +  IWIDTH  FFT00800 

TMP  =  X(IBOT+l)*WTFAC  FFT00810 

XCIBOT+l)  =  XCITOP+1)  -  TMP  FFT00820 

X( ITOP+1)  =  X(ITOP+l)  +  TMP  FFT00830 

30  CONTINUE  FFT00840 

40  CONTINUE  FFT00850 

50  CONTINUE  FFT00860 

RETURN  FFT01000 

END  FFT01010 
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5.  \RDy\* 

c** 

THIS  FORTRAN  FILE  COMPUTES  THE  SQUARED  MAGNITUDE  OF  THE  *** 

C** 

RIHACZEK  DISTRIBUTION 

ititit 

c 

* 

c 

* 

♦ 

c 

INPUT  DATA  SEQUENCE  IS  READ  USING  FILEDEF  4,  AS  THE  COMPLEX 

* 

c 

ARRAY  X(N). 

* 

c 

* 

c 

N  IS  THE  LENGT.i  OF  THE  DATA  SEQUENCE  AND  IS  ADJUSTED  FROM  THE 

* 

4 

c 

PARAMETER  STATEMENT. 

N  MUST  NOT  EXCEED  128. 

* 

c 

* 

c 

ANALYSIS  PARAMETERS  ARE  READ  USING  FILEDEF  41.  THE  PARAMETERS 

* 

c 

ARE: 

* 

c 

ARGUMENT  TYPE 

ALLOWED  VAULUES 

* 

c 

* 

c 

MODE  -  11 

1  PLOT  0  TO  PI 

* 

c 

2  PLOT  PI  TO  PI 

* 

\ 

c 

it 

c 

PLTR  -  11 

0  SHERPA  LASER  PRINTER 

* 

l 

c 

1  IMB79  GRAPHICS  TERMINAL 

* 

i 

c 

* 

\ 

c 

BWLEN-  13 

3  DIGIT  INITIAL  WINDOW  LENGTH, 

* 

< 

c 

MUST  BE  AN  ODD  INTEGER 

* 

c 

EWLEN-  13 

3  DIGIT  FINAL  WINDOW  LENGTH 

* 

c 

WINC  -  12 

2  DIGIT  WINDOW  INCREMENT,  MUST 

* 

1 

c 

BE  AN  EVEN  INTEGER 

* 

c 

* 

•  '■ 

c 

WTYPE-  A19 

19  CHARACTER  STRING  USED  IN  THE 

* 

c 

PLOT  HEADER  DISCRIBING  THE 

* 

c 

WINDOW  USED.  THE  CURRENT 

* 

1 

c 

WINDOW  LENGTH  IS  AUTOMATICALLY 

it 

c 

INCLUDED 

it 

c 

* 

c 

CONTR-  A1 

1  CHARACTER  STRING  INDICATING 

* 

c 

TYPE  OF  PLOT  DESIRED 

it 

c 

it 

c 

A  AMPLITUDE  PLOT  ONLY 

* 

c 

C  CONTOUR  PLOT  ONLY 

* 

c 

B  BOTH  AMPLITUDE  AND  CONTOUR 

* 

■; 

c 

* 

5 

c 

TTL  -  A43 

43  CHARACTER  STRING  USED  IN  THE 

* 

j 

c 

HEADING  WHICH  DESCRIBES  THE 

* 

1 

1  c 

ALGORITHM  AND  THE  CLASS  OF 

* 

c 

SIGNAL  USED 

* 

c 

it 

c 

SIGNAL-  A43 

43  CHARACTER  STRING  DESCRIBING 

it 

c 

TEST  SIGNAL 

it 

\ 

c 

it 

\ 

c 

TSMTH-  A25 

25  CHARACTER  STRING  DESCRIBING 

it 

\ 

c 

TYPE  OF  TIME  SMOOTHING  USED 

* 

3, 

i 

c 

it 

1 

.  ) 

c 

* 

j 

c 

OUTR  REAL 

it 

c 

INITIALLY  USED  AS  STORAGE  FOR  THE  TIME  SMOOTHED  IPS,  THEN 

it 

c 

AS  OUTPUT  ARRAY  OF 

DIMENSION  512  BY  N 

it 
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* 

OUTI  REAL  * 

ARRAY  USED  TO  HOLD  THE  THE  TIME  SMOOTHED  IMRD  * 

* 

SAMP  COMPLEX  * 

SHIFTED  VERSION  OF  X  * 

* 

SAM  COMPLEX  * 

SHIFTED  AND  CONJUGATED  VERSION  OF  X  * 

* 

RE  COMPLEX  * 

ARRAY  OF  SUM  OF  PRODUCTS  OF  DIMENSION  1  BY  1025,  * 

POSITIONS  1  -  512  ARE  CONJUGATE  SYMMETRIC  WITH  * 

POSITIONS  514  -  1025,  POSITION  513  =  (0. ,0. )•  * 

* 

IM  COMPLEX  * 

ARRAY  OF  DIFFERENCE  OF  PRODUCTS  OF  DIMENSION  1  BY  1025,  * 

POSITIONS  1  -  512  ARE  CONJUGATE  SYMMETRIC  WITH  * 

POSITIONS  514  -  1025,  POSITION  513  =  (0. ,0. )•  * 

* 

FTR  COMPLEX  * 

ARRAY  OF  THE  1024  POINT  TRANSFORM  RE  * 

* 

FTI  COMPLEX  * 

ARRAY  OF  THE  1024  POINT  TRANSFORM  IM 

* 

Z  INTEGER  * 

LENGTH  OF  CURRENT  WINDOW  * 

* 

M  INTEGER  * 

MID -POINT  OF  THE  CURRENT  WINDOW  * 

* 

AMAX  REAL  * 

MAXIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 

AMIN  REAL  * 

'  MINIMUM  AMPLITUDE,  USED  TO  SCALE  VERTICAL  AXIS  * 

* 

NX  HORIZONTAL  DIMENSION  OF  OUT  WHICH  IS  ALWAYS  512  * 

* 


DATA  CAN  BE  OUTPUT  USING  FILEDEF  61.  POINTS  OF  INTEREST  MUST  * 
BE  DEFINED  IN  THE  APPROPRAITE  SECTION  OF  CODE.  BECAUSE  * 
OF  SPACE  CONSTRAINTS,  THE  DATA  OUTPUT  FILE  IS  WRITTEN  TO  * 
THE  B  DISK.  * 


PARAMETERS  64) 

COMPLEX  X(N) ,RE( 1025) , SAMP, SAM, FTR( 1024) 

COMPLEX  IM( 1025), FTI( 1024) 

REAL  OUTR(512,N),OUTI(512,N), AMAX, AMIN 
INTEGER  NX , K , 1 , J , MODE , Z , M , BWLEN . EWLEN , WINC , PLTR 
CHARACTER  WTYPE*19 ,TTL*43 , SIGNAL*37 ,TSM7H*25 , C0NTR*1 
CALL  EXCMSC FILEDEF  4  DISK  TEST  IN  (PERM') 

CALL  EXCMS( 'FILEDEF  41  DISK  PARAM  IN  (PERK  ) 
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CALL  EXCMS( 'FILEDEF  61  DISK  DATA  OUT  B  (PERM’) 

■ . READ  IN  PARAMETER  LISITNG . - . 

READ  (41,400) MODE , PLTR , BWLEN , EWLEN , WINC , WTYPE , CONTR , TTL , 

+  SIGNAL, TSMTH 

FORMAT  ( 1X,I1,3X,I1,3X,I3,3X,I3,3X,I2/1X,A19/1X,A1/1X,A43 
+  /1X,A37/1X,A25) 


.  TEST  TO  ENSURE  WINDOW  LENGTH  IS  APPROPRIATE  . 

I=N-1 

IF  ((BWLEN  .GT.  I)  .OR.  (EWLEN  .  GT.  I))  THEN 
WRITER, 69) 

GO  TO  99 
ENDIF 

I=MOD(BWLEN,2) 

K=MOD(WINC,2) 

IF  (I  .EQ.  0)  THEN 

IF  (K  .EQ.  1)  THEN 
WRITE(*,68) 

GO  TO  99 

ELSE 

WRITE(*,67) 

GO  TO  99 
ENDIF 

ENDIF 

FORMAT  (IX,’ WINDOW  LENGTH  EXCEEDS  LENGTH  OF  THE  DATA’) 
FORMAT  (IX, ’WINDOW  INCREMENT  MUST  BE  EVEN’) 

FORMAT  (IX, 'INITIAL  WINDOW  LENGTH  MUST  BE  ODD') 


.  PLOTTING  DEVICE  CALL 

IF  (PLTP.EQ.O)THEN 
CALL  COMPRS 

ELSE 

CALL  TBM79 
ENDIF 


PI=4*ATAN( 1.  ) 

NX=512 

READ  (4,*)(X(J),J=1,N) 

DO  111  Z=  BWLEN, EWLEN, WINC 
WRITE  (61, 600)rrL, SIGNAL, TSMTH, WTYPE, Z 
600  FORMAT  ( 1X,A43/1X,A37/1X,A25/1X,A19,I3, '  POINTS)') 
CALL  ANGLE (0.  ,0.  ) 

M=(Z-l)/2 

AMAX=0. 

AM I N= AM AX 
DO  10  1=1, N 
DO  20  K=0,512 
SAMP=(0.  ,0. ) 

SAM=SAMP 

IF  (  (I+K)  .LE.  N  )  THEN 
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SAMP=X( I+K) 

ENDIF 

IF  (  (I-K)  .GT.  0  )  THEN 
SAM=CONJG( X( I -K) ) 

ENDIF 

C 

IF  (K  .IE.  M)  THEN 

RE(K+1)=(X(I)*SAM  +  CONJG( X( I ) )*SAMP ) 

+  *( 0. 54+0.  46*COS( 2*PI*K/( 2*M) ) ) 

C 

IM(K+1)=(X(I)*SAM  -  CONJG( X( I ) )*SAMP) 

+  *( 0. 54+0. 46*COS( 2*PI*K/( 2*M) ) ) 

ELSE 

RE(K+1)=0 

IM(K+1)=0 

ENDIF 

C 

RE(1024-K+1)=CONJG(RE(K+1)) 

IM( 1024-K+1)=C0NJG( IM(K+1) ) 

20  CONTINUE 

RE(513)=(0.  ,0.  ) 

IM(513)=(0.  ,0.  ) 

CALL  FFT(1024,RE,FTR) 

CALL  FFT( 1024,IM,FTI) 

DO  40  K=l, M0DE*5 12, MODE 
C 

IF  (  (K  .LT.  514)  .AND.  (MODE  .EQ.  2))  THEN 
OUTR( INT((K+l)/2+255) ,I)=REAL(FTR(K)) 

OUTI ( INT( ( K+l ) / 2+255 ) , I )=REAL( FTI( K) ) 

ELSE 

IF  (  MODE  .EQ.  2  )  THEN 

OUTR( INT( ( K+ 1 ) / 2 -25  7 ) , I ) -REAL( FTR(K)) 
OUTI(INT((K+l)/2-257) ,I)=REAL(FTI(K)) 

ELSE 

OUTR( K , I )=REAL(FTR(K) ) 

OUTI ( K , I )=REAL( FTI ( K) ) 

ENDIF 

ENDIF 

40  CONTINUE 

10  CONTINUE 

C 

C . FOR  TIME  SMOOTHING . 

DO  48  K=l,512 
DO  46  1=1, N-2 

OUTR( K , I )=( OUTR( K , I ) +OUTR( K , 1+1 )+OUTR( K , 1+2 ) ) / 3 
OUTI(K,I)=(OUTI(K,I)+OUTI(K,I+l)+OUTI(K,I+2))/3 

46  CONTINUE 

DO  47  I=N,3,-1 

OUTR(K, I)=(OUTR(K, I)+OUTR(K, I-l)+OUTR(K, 1-2) )/3 
OUTI(K,  I)=(  OUTI(K,  I')+OUTI(K,  I-l)+OUTI(K,I-2)  )/3 

47  CONTINUE 

48  CONTINUE 


C+  -+-+“+-+-+-+“+-+-+"  -+-+-+-+-+-+-+ 

C+-+-  the  sura  of  magnitudes  is  |RD|**2  Difference  is  |RD-|**2  +-+ 
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DO  200  K=1 .512 
DO  201  1=1, N 

OUTR( K , I )=ABS( OUTR( K , I ) ) 

OUTI(K,I)=ABS(OUTI(K,I)) 

OUTR(K,I)=OUTR(K,I)  +  OUTI(K,I) 

C 

IF  (OUTR(K,I)  .GT.  AMAX)  THEN 
AMAX=OUTR(K,I) 

ENDIF 

IF  (OUTR(K,I)  .LT.  AMIN)  THEN 
AMIN=OUTR(K,I) 

ENDIF 

201  CONTINUE 
200  CONTINUE 

C+  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  .  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  + 


50 


PLOTTING  . 

IF  (CONTR.EQ.  ’C’)THEN 
GO  TO  50 
ENDIF 

CALL  HEIGHT(0.  28) 

CALL  BSHIFT  (  -0.2,-* 25) 

CALL  AREA2D(8,9) 

CALL  VOLM3D(10,10,8) 

CALL  HEADIN(TTL, 100,1.  ,3) 

CALL  HEADIN( SIGNAL, 100,1. ,3) 

CALL  HEADIN(TSMTH, 100,1. ,3) 

CALL  ME3SAG(WTYPE, 100.2.  5,9.  3) 

CALL  INTNOCZ  , ' ABUT' , 'ABUT' ) 

CALL  MESSAGE  POINTS)  $’  ,100/  ABUT' ,’  ABUT') 

CALL  X3NAME( 'FREQUENCY  AXIS$',100) 

CALL  Y3NAME( 'TIME  AXIS?', 100) 

CALL  Z3NAME( '  $',100) 

CALL  VUANGL( -65 ,70,700) 

CALL  XNONUM 
CALL  ZNONUM 

CALL  MX1ALF( ’ STANDARD ','#') 

CALL  MX2ALF( ' L/CGREEK' , ’+' ) 

CALL  ANGLE (-25.  0) 

IF  (  MODE  .EQ.  2  )  THEN 
CALLMESSAGC’  +-P#  ’,6,0. ,2.3) 

ELSE 

CALL  MESSAGC  +0#  ',5,0.  ,2.3) 

ENDIF 

CALL  ANGLE (-25.  0) 

CALLMESSAGC  +P#  ’,5,4.9,0.15) 

CALL  GRAF3D( -256,256,256, 1,N,N, 1. 0*AMIN,0.  5*(AMAX-AMIN) , 
+  1. 0*AMAX) 

CALL  SURMAT(0UTR,512,512, 1,N,0.  ) 

CALL  ENDPL(O) 

CONTINUE 

IF  (CONTR.NE. 'A')THEN 

CALL  CONTOR( OUTR , NX , N , TTL , SIGNAL , WTYPE , TSMTH , Z ) 

ENDIF 
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WRITE(61,603)AMAX,AMIN 
603  FORMAT  ( IX,  ’ MAXIMUM  AMPLITUDE*2’  ,E14.  7 
+  /IX, 'MINIMUM  AMPLITUDE2*'  ,E14.  7) 

111  CONTINUE 
CALL  DONEPL 
99  STOP 
END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 

SUBROUTINE  C0NT0R( A , NX ,NY , TITLE , S IGNL , VrNDW ,TAVG , WLEN) 

THIS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS. 
NOTE:  THE  ARRAY  MUST  BE  REAL*4. 

A  :  SINGLE  PRECISION  NX  BY  NY  ARRAY  OF  REGULARLY  SPACED  POINTS 
NX:  NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NY:  NUMBER  OF  POINTS  IN  THE  Y -DIRECTION 
ZINC:  CONTOUR  INTERVAL 

DIMENSION  A(NX,NY) 

COMMON  W0RK(50000) 

SET  PARAMETERS  FOR  AXES: 

XORIG=-256. 

XSTP=256. 

XMAX=256. 

YORIG=0. 

YSTP=NY 

YMAX=NY 

SET  CONTOUR  LEVEL 
ZINC=AMAX/10 

CALL  SETCLR('CYAN') 

SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT: 

CALL  PAGE( 8.  5,11.0) 

CALL  BC0M0N( 50000) 

CALL  AREA2D(6.  0,8.  0) 

LABEL  AXES: 

CALL  XNAMEC 'FREQUENCY  -  AXIS  $',100) 

CALL  YNAME( ' TIME  -  AXIS  $',100) 

CALL  XNONUM 

CALL  GRAF (XORIG, XSTP , XMAX , YORI G , YSTP , YMAX ) 

CALL  FRAME 

TITLE: 

CALL  HEADINC ' CONTOUR  PLOT? ', 100, 1.  ,4) 

CALL  HEADIN(TITLE, 100,1. ,4) 

CALL  HEADINC SIGNL, 100,1.  ,4) 

CALL  HEADIN(TAVG, 100,1.  ,4) 

CALL  ANGLE (0.  0)  ■  ’ 

CALL  MESSAG(WNDW, 100,1. 5,-.  7) 

CALL  INTNOCWLEN  , 'ABUT' , 'ABUT' ) 
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CALL  MESSAG( ' POINTS) $ ' , 100 , ' ABUT' , ' ABUT' ) 

C 

C  MAKE  CONTOURS  AND  DRAW: 

CALL  SETCLR('RED') 

CALL  C0NMIN(3. 0) 

CALL  CONANG(60.  ) 

CALL  C0NLIN(0,'MYCON' ,’ NOLABELS' ,2,10) 

CALL  C0NMAK( A , NX , NY , ’ SCALE ' ) 

CALL  C0NTUR( 1 , ' LABELS ’ , ' DRAW ' ) 

C 

C  END  PLOT: 

CALL  ENDPL(O) 

RETURN 
END 

ssssssssssssssssssssssssssssssssss 

SUBROUTINE  MYCON(RARAY, IARAY) 

DIMENSION  RARAY(2),IARAY(1) 

THIS  ROUTINE  MAKES  NEGATIVE  CONTOURS  DASHED  AND  THE  ZERO  LINE  HEAVIER. 

CALL  RESET( 'DASH' ) 

IF  (RARAY(l)  .GE.  0.)  GO  TO  10 
CALL  DASH 
10  RARAY(2)  =  1. 

IARAY(l)  =  1 

IF  (RARAY(l)  .EQ.  0.)  IARAY(l)  =  2 
RETURN 
END 

SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS 


* 


* 


*  CALL  FFT(N,XTEMP,X) 

* 


* 

* 


*  X  -  OUTPUT  COMPLEX  ARRAY  CONTAINING  FFT  (1024)  * 

*  N  -  NUMBER  OF  POINTS  * 

*  XTMP  -  COMPLEX  ARRAY  CONTAINING  DATA  SAMPLES  * 

*  (starting  at  l,up  to  1024)  * 

******yf**********Vf***yf*********!-ryf******yr***********Vf*Vf 


SUBROUTINE  FFT(N,XTMP,X) 

COMPLEX  X( 1024) ,XTMP( 1024) ,WTFAC,TMP 
M=INT(LOGlO(FLOAT(N))/LOG10(2. )+0. 5) 

EN  =  N 

PI  =  4. 0*ATAN(1. 0) 

DO  10  K=0,N-1 
NEWADR  =  0 
MADDR  =  K 
DO  20  1=0, M-l 

LRMNDR  =  MOD (MADDR, 2) 

NEWADR  =  NEWADR  +  LRMNDR*2**(M-1-I) 
MADDR  =  MADDR/ 2 
20  CONTINUE 


FFT00130 

FFT00140 

FFT00210 

FFT00270 

FFT00320 

FFT00330 

FFT00340 

FFT00350 

FFT00360 

FFT00370 

FFT00380 

FFT00390 
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X(NEWADR+1)  =  XTMP(K+1) 

10  CONTINUE 
DO  50  L=1 ,M 

ISPACE  =  2**L 
S  =  N/ ISPACE 
IWIDTH  =  I SPACE/ 2 
DO  40  J=0,( IWIDTH- 1) 

R  ss 

ALPHA  =  2.*PI*R/EN 

WTFAC  =  CMPLX(  COS( ALPHA),  -SIN( ALPHA)) 
DO  30  IT0P=J,N-2, ISPACE 
IBOT  =  ITOP  +  IWIDTH 
TMP  =  X(IBOT+l)*WTFAC 
X( IBOT+1)  =  X(ITOP+l)  -  TMP 
XdTOP+l)  =  X(IT0P+1)  +  TMP 
30  CONTINUE 

40  CONTINUE 
50  CONTINUE 
RETURN 
END 


FFT00400 

FFT00410 

FFT00530 

FFT00610 

FFT00620 

FFT00630 

FFT00670 

FFT00720 

FFT00730 

FFT00740 

FFT00750 

FFT00800 

FFT00810 

FFT00820 

FFT00830 

FFT00840 

FFT00850 

FFT00860 

FFT01000 

FFT01010 
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APPENDIX  B.  CROSS  IPS 


Up  to  this  point  only  autospectra  have  been  discussed.  Analysis  of  cross  spectral 
characteristics  of  nonstationary  phenomena  can  provide  valuable  information  about  the 
process.  Equivalent  to  (1),  a  cross  power  spectral  density  is  defined 


Pv,V)-r  Rv,Wn’fi<*r- 

— oo 


(63) 


For  the  case  when  x,  and  Xj  are  uncorrelated  then 

Px^jJ)  -  2 nninj5{f) 


(64) 


where  n,  is  the  mean  value.  If  the  data  is  correlated  the  energy  resulting  from  cross 
spectral  analysis  can  be  complex.  By  examining  Parseval's  theorem  in  a  more  general 
context, 


=i v  (65) 

where  R  ( 0)  is  not  necessarily  real  nor  is  the  cross  correlation  function  (CCF)  neces¬ 
sarily  conjugate  symmetric  about  R  ( 0).  [Ref.  1] 

All  the  spectral  estimators  previously  discussed  are  applicable  if  the  ACF  estimate 
is  replaced  by  with  the  CCF  estimate.  The  bias  for  cross  spectra  may  be  much  larger 
than  an  equivalent  autospectra  where  the  point  of  maximum  overlap  occurs  at  lag  zero. 
In  practice,  the  location  of  maximum  overlap  is  unknown  for  CCF's.  One  interpretation 
of  cross  PSD  are  of  importance,  the  case  where  x,  and  Xj  are  two  channels  of  a  multi¬ 
channel  system.  The  \PvJf)  I  contains  information  concerning  relative  amplitudes  at 
specific  frequencies  where  the  4^  contains  information  concerning  the  lead  or  lag  in 
phase  between  the  two  channels.  [Ref.  3] 
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Applying  IPS  in  a  multi-sensor  environment  leads  to  the  following  defining 
equation, 


1*00 


/PS^,0  =  t 


(*(*)/(*  -  t)  +  x\t)y(t  +  t))  e~Jlnj‘dxt 


(66) 


where  signal  y  could  be  a  delayed,  noisy  version  of  signal  x.  Using  a  rapidly  changing, 
linearly  chirped  pulse  and  implementing  a  windowed  version  of  IPS,#, 

two  sample  cross  spectra  are  shown  in  Figure  50  and  Figure  51. 

The  first  spectra  was  created  by  beating  the  pulsed  chirp  x{t )  against  a  delayed  ver¬ 
sion  of  itself.  The  cross  spectrum  is  shown  in  Figure  50  (a).  In  this  case  the  delay  is 
18  samples  for  a  pulse  63  samples  in  duration.  The  cross  spectrum  can  be  seen  in 
Figure  50  (b)  and  (c),  where  the  spectral  ridge  over  the  interval  corresponding  to  the 
absolute  time  of  overlap.  The  maximum  amplitude  achieved  on  the  cross-spectral  sur¬ 
face  is  nearly  the  same  as  for  the  autospectra.  The  minimum  however,  is  approximately 
44%  greater  in  magnitude  than  that  found  on  the  corresponding  auto-spectral  surface. 
1PSXj  does  not  appear  to  provide  information  which  can  estimate  delay  in  reception  for 
this  class  of  signal. 

The  second  cross  spectrum  considered  examines  the  ability  of  IPSXJt  to  indicate  cor¬ 
relation  between  a  pulsed  chirp  and  a  Doppler-shifted  version  of  itself.  In  this  case, 

.*(/)  =  e^K(  lW f+20(  T2F )J) 

/is  t  2\  (67) 

y(t)  =  128  f+20(  128  )  ), 


representing  a  shift  of  50%.  Figure  51  (c)  shows  an  overlay  of  the  two  pulses. 
Figure  51  (b)  is  the  contour  for  the  cross  spectrum.  The  peaks  of  the  characteristically 
modulated  ridge  correspond  to  the  line  of  overlap  shown  in  (c).  It  is  not  clear  if  the  cross 
spectrum  of  a  linear  chirp  with  a  Doppler  shifted  version  of  itself  yields  information 
concerning  the  degree  of  coherence.  The  cross  spectrum  in  Figure  51  could  easily  be 
interpreted  as  an  autospectrum  in  which  two,  closely -spaced  parallel  chirps  are  present. 
This  initial  investigation  into  the  behavior  of  IPS TOf  suggests  that  a  more  detailed  exam¬ 
ination  of  its  behavior  is  in  order. 
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(c)  IPSXJ  amplitude  plot 

Figure  50.  Cross  spectral  analysis  of  a  pulsed  linear  chirp 
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(c)  IPS,  overlayed  contour  plots  of  the  original  and  Doppler-shifted  pulse 


Figure  51.  Cross  spectral  analysis  for  a  Doppler-shifted  linear  pulse 


APPENDIX  C.  CUMULANT 


Assuming  the  received  signal  is  corrupted  by  zero  mean,  Gaussian  noise,  examining 
the  third-order  moment  or  cumulant  may  yield  information  about  the  signal  while  sup¬ 
pressing  the  contributions  of  the  noise.  This  potential  processing  gain  is  realized  because 
the  odd  moments  of  a  Gaussian  process  are  identically  zero.  One  way  to  implement  IPS 
as  an  estimator  of  the  cumulant  is 


/•  oo 

Cum.  IPSf  KD  |  (f,t)  -  y  ( I  *(/)  1 2jr(r  - 1)  +  I  *'(/)  I  !jt(i  +  t)) 


Looking  only  at  the  effect  on  noiseless  signal  data,  some  preliminary  results  can  be 
seen  in  Figure  52.  All  signals  are  analytic,  therefor  the  magnitude  squared  term  is  al¬ 
ways  unity.  A  comparison  of  the  treatment  of  IPSy  can  be  made  by  referring  to  Chapter 
IV,  Section  C  (3):  Test  Case  Results.  This  method  of  processing  the  cumulant  of  an 
analytic  signal  does  not  appear  to  provide  any  useful  results.  Further  research  should 
examine  the  behavior  using  real  signals  or  possibly  forming  the  approximation  of  the 
cumulant  in  a  different  fashion. 


(d)  FSK 


Figure  52.  The  cumulant  of  various  analytic  signals 


107 


TUC  «1S 


LIST  OF  REFERENCES 


1.  P.  B.  Peebles,  Jr.,  " Spectral  characteristics  of  random  processes,"  in  Probability, 
Random  Variables  and  Random  Signal  Principles ,  pp.  172-204,  McGraw  Hill,  New 
York,  1980. 

2.  Athanasios  Papoulis,  " Correlation  and  power  spectrum  of  stationary  processes,"  in 
Probability,  Random  Variables  and  Stochastic  Processes  ,  pp.  336-384,  McGraw 
Hill,  New  York,  1965. 

3.  Steven  M.  Kay,  Modern  Spectral  Estimation:  Theory  and  Application , 

McGraw  Hill,  New  York,  1980  • 

4.  Fredric  J.  Harris,  "On  the  Use  of  Windows  for  Harmonic  Analysis  with  the  Discrete 
Fourier  Transform,"  IEEE  Proc.,  Vol.  66,  No.  1  pp.  51-83,  January  1978. 

5.  M.  H.  Ackroyd,  "Instantaneous  and  Time-Varying  Spectra:  an  Introduction,"  The 
Radio  and  Electronic  Engine  ,  Vol.  39,  No.  3,  pp.  142-152,  March  1980. 

6.  O.  D.  Grace,  "Instantaneous  Power  Spectra,"  J.  Acoust.  Soc.  Amer.,  Vol.  39,  No. 
3,  pp.  191-197,  January  1981. 

7.  W.  Mecklenbrauker,  "A  Tutorial  on  Non-Parametric  Bilinear  Time- Frequency  Sig¬ 
nal  Representations,"  Les  Houches,  Session  XLV,  1985,  J.L.  Lacoume  and  R.  Stora, 
eds.,  pp.  276-236,  Elsevier  Science  Publishers  B.V.,  1987. 

8.  Richard  B.  Altes,  "Detection,  Estimation,  and  Classification  with  Spectrograms," 
J,  Acoust,  Soc.  Amer.,  Vol.  67,  No.  4,  pp.  1232-1246.  April  1980. 

9.  Leon  Cohen,  "Generalized  Phase-Space  Distribution  Functions,"  J.  Math  Physics, 
Vol.  7,  No.  5,  pp.  781-786,  May  1966. 


108 


% 


< 


10.  Leon  Cohen,  "Time-frequency  Distributions  -  a  Review,"  Proc .  IEEE,  Vol.  77,  No. 
7,  pp.  941-981,  July  1989- 

11.  C.  H.  Page,  "Instantaneous  Power  Spectra,"  J  Appl.  Phys.,  Vol.  23,  pp.  103-106, 
1952. 

12.  Morris  J.  Levin,  "Instantaneous  Spectra  and  Ambiguity  Functions,"  Trans.  IEEE 
on  Information  Theory ,  Vol.  IT-10,  pp.  95-96,  January  1 964  • 

13.  R.M.  Fano,  "Short-Time  Autocorrelation  Functions  and  Power  Spectra,"  J.  Acoust. 
Soc.  Amer.,  Vol.  22,  No.  5,  pp.  546-550,  September  1950. 

14.  M.R.  Schroeder  and  B.S.  Atal,  "Generalized  Short-Time  Power  Spectra  and 
Autocorrelation  Functions,"  J.  Acoust.  Soc.  Amer.,  Vol.  34,  No.  11,  pp.  1679-1683, 
November  1962. 

15.  Paulo  M.  D.  Monica  de  Oliveira,  "Instantaneous  Power  Spectrum Engineer's  The¬ 
sis,  Naval  Postgraduate  School,  Monterey,  California,  March  1989. 

16.  August  W.  Rihaczek,  "Signal  Energy  Distribution  in  Time,"  Trans.  IEEE  on  Infor¬ 
mation  Theory,  Vol.  14,  No.  3,  pp.  369-374,  May  1968. 

17.  E.  P.  Wigner,  "On  the  Quantum  Correction  for  Thermodynamic  Equilibrium,"  Phys. 
Rev.,  Vol.  40,  No.  pp.  749-759,  1932. 

18.  J.  Ville,  "Theorie  et  Applications  de  la  Notion  de  Signal  Analytique,"  Cables  et 
Transmission,  Vol.  2,  No.  pp.  61-74,  1948. 

19.  T.A.C.M.  Claasen  and  W.F.G.  Mecklenbrauker,  "  The  Wigner  Distribution  -  a 
Too!  for  Time  Frequency  Signal  Analysis,  Part  III:  Relations  with  other  Time- 
Frequency  Signal  Transforms,"  Phillips  J.  Research,  Vol.  35,  No.  6,  pp.  372-389, 
1980. 


109 


20.  Boualem  Bouashash  "Notes  on  the  use  of  the  Wigner  Distribution  for  Time- 
Frequency  Signal  Analysis,"  Trans.  IEEE  ASSP ,  Vol.  36,  No.  9,  pp.  1518-1521, 
September  1988- 

21.  Cornelis  P.  Janse  and  Arie  J.M.  Kaiser  "Time-Frequency  Distribution  of  Loud¬ 
speakers:  the  Application  of  the  Wigner  Distribution,"  J.  Audio  Eng.  Soc .,  Vol.  31, 
No.  4,  pp.  198-216,  April  1983- 

22.  H.I.  Choi  and  W.J.  Williams,  "Improved  Time-Frequency  Representations  of 
Multicomponent  Signals  Using  Exponential  Kernels,"  Trans  IEEE  on  Acoust., 
Speech  and  Signal  Processing,  Vol.  ASSP-37,  1989* 


110 


BIBLIOGRAPHY 


Allard,  Jean  F.  and  Valiere,  Jean  C.  and  Bourdier  R.,  "Broadband  Signal  Analysis 
with  the  Smoothed  Pseudo-Wigner  Distribution,"  J.  Acoust.  Amer.,  Vol.  83,  No.  3, 
pp.  1041-1044,  March  1988. 

Amin,  Moeness  G.,  "Sliding  Spectra:  A  New  Perspective,"  Fourth  Annual  ASSP 
Workshop  on  Spectrum  Estimation  and  Modeling ,  pp.  55-59,  IEEE  Press,  New  York, 
1988. 


T.A.C.M.  Claasen  and  W.F.G.  Mecklenbrauker,  "The  Wigner  Distribution  -  a  Tool 
for  Time  Frequency  Signal  Analysis,  Part  I:  Continuous  Time  Signals,"  Phillips  J . 
Research,  Vol.  35,  No.  6,  pp.  217-250,  1980. 

T.A.C.M.  Claasen  and  W.F.G.  Mecklenbrauker,  "The  Wigner  Distribution  -  a  Tool 
for  Time  Frequency  Signal  Analysis,  Part  II:  Discrete  Time  Signals,"  Phillips  J. 
Research,  Vol.  35,  No.  6,  pp.  276-300,  1980. 

Cohen,  Leon  and  Posch,  Theodore  E.,  "Positive  Time-Frequency  Distribution 
Functions,"  Trans.  IEEE  on  Accoustics,  Speech  and  Signal  Processing,  Vol.  ASSP 
33,  No.  1,  pp.  31-38,  February  1985. 

Flaska,  M.,D.,  "Cross  Correlation  of  Short  Time  Spectral  Histories,"  Acoust.  Soc. 
Amer. ,  Vol.  59,  No.  2,  pp.  381-388,  February  1976. 

Lampard,  D.G.,  "Generalization  of  the  Wiener- Khintchine  Theorem  to  Nonsta¬ 
tionary  Processes,"  J.  App.  Physics,  Vol.  25,  No.  6,  pp.  802-803,  June  1954. 


Mark,  W.,D.,  "Spectral  Analysis  of  the  Convolution  and  Filtering  of  Nonstatio 
Stochastic  Processes,"  J.  Sound  Vib.,  Vol.  11,  pp.  19-63,  1970. 


nary 


111 


Martin,  W.  and  Flaodrin,  P.,  "Analysis  of  Nonstationary  Processes  Short  time 
Periodograms  versus  a  Pseudo-Wigner  Estimator,"  Signal  Processing  11:  Theories 
and  Applications,  edited  by  H.W.  Schtlssler,  pp.  455-458,  Elsevier  Science  Publishers 
B.V.,  1983. 

Stutt,  Charles  A.,  "Some  Results  on  Real-Part/Imaginary-Part  and 
Magnitude/Phase  Relations  in  Ambiguity  Functions,"  Trans.  IEEE  on  Information 
Theory ,  pp.  321-327,  October  1964. 

Urkowitz,  Harry,  "Pre-envelopes  of  Nonstationary  Bandpass  Processes,"  J.  Franklin 
Institute,  Vol.  277,  No.  1,  pp.  31-36,  January  1964. 


112 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 


1.  Defense  Technical  Information  Center  2 

Cameron  Station 

Alexandria,  VA  22304-6145 

2.  Library,  Code  0142  2 

Naval  Postgraduate  School 

Monterey,  CA  93943-5002 

3.  Professor  Ralph  D.  Hippenstiel,  Cod?  EC/Hi  3 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA,  93943-5004 

4.  Professor  Roberto  Cristi.  Code  EC/Cx  1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postsraduate  School 
Monterey,  CA,  93943-5004 

5.  Chairman,  Code  EC  1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA,  93943-5004 

6.  Professor  Charles  W.  Therrien,  Code  EC/Th  1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA,  93943-5004 

7.  Professor  Murali  Tummala,  Code  EC/Tu  1 

Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA,  93943-5004 

8.  Naval  Ocean  Systems  Center  1 

Attn:  Dr.  C.  E.  Persons,  Code  732 

San  Diego,  CA,  92152-5000 

9.  Naval  Ocean  Systems  Center  3 

Attn:  LT  E.K  Stitz,  Code  30 

San  Diego,  CA,  92152-5000 


113 


